PDF
Solid Oxide Electrolyzer Using Thermodynamics
Introduction
This example models a solid oxide electrolyzer cell wherein water vapor is reduced to form hydrogen gas on the cathode, and oxygen gas is evolved on the anode. The current distribution in the cell is coupled to the cathode mass transfer of hydrogen and water and momentum transport.
Model Definition
On the anode, oxygen ions are oxidized to form oxygen gas,
(1)
whereas on the cathode, water vapor is reduced to form hydrogen gas and oxygen ions:
(2)
Figure 1 shows the model geometry. Since the oxygen is the only gas present in the anode gas chamber, and isobaric conditions are assumed, there is no need to explicitly model the anode gas transport. Four computational domains are hence used in the model: the cathode gas channels, the cathode gas diffusion electrode, the solid oxide electrolyte layer, and the anode gas diffusion electrode.
Figure 1: Model geometry. From top: Cathode gas channels, cathode gas diffusion electrode, solid oxide electrolyte layer, and anode gas diffusion layer. The positions of the inlet and outlet are indicated in the figure.
The composition of the hydrogen-water vapor mixture will change as a result of the electrochemical reactions. The mass transport of hydrogen and water vapor is modeled in the cathode gas channels and the gas diffusion electrode, coupled to the resulting (laminar) flow of the gas mixture. The mass transport of hydrogen and water is modeled using the Transport of Concentrated Species interface, using Maxwell-Stefan diffusion. The momentum flow is defined in the model using the Brinkman Equations interface for the porous gas diffusion electrodes. The Navier-Stokes equations are used for the nonporous gas channels.
The current distribution is defined assuming a constant conductivity of the solid electrolyte. The Secondary Current Distribution interface is used to define the electrode reactions and the electrolyte charge transport in the porous gas diffusion electrodes and the electrolyte layer.
On the cathode side, the electrode kinetics depends on the local concentration of water and hydrogen according to the law of mass action (and Nernst equation). On the anode side, and a uniform partial pressure of oxygen is assumed and a concentration-independent Butler-Volmer expression is hence used to define the electrode kinetics.
The Thermodynamics and Chemistry nodes are used to automatically define the properties of the cathode gas mixture, as well as the equilibrium potentials of the electrode reactions.
Results and Discussion
Figure 2 shows the velocity magnitude distribution in the cell. The highest velocities are located close to the inlet and outlet.
Figure 2: Velocity in the cell.
Figure 3 shows how the density and the dynamic viscosity of the gas relate to the hydrogen and water molar fractions shown in Figure 4. As the hydrogen content of the gas increases toward the outlet, the density and the viscosity both decrease.
Figure 3: Density (left) and dynamic viscosity (right).
Figure 4: Hydrogen (left) and water vapor (right) concentrations.
Figure 5 shows the molar fraction of hydrogen in the gas mixture, and the corresponding hydrogen flux streamlines. The molar fraction is close to zero at the inlet and surpasses 80% at the outlet.
Figure 5: Hydrogen molar fraction (slice) and flux (streamlines).
Figure 6 shows the cross-sectional electrolyte current density in the middle of the electrolyte between the anode and the cathode. The current density is highest close to the inlet, where the water/hydrogen ratio is high, and decreases toward the outlet.
Figure 6: Electrolyte current density through electrolyte layer.
Finally, Figure 7 shows the cell voltage as a function of the average current density (polarization curve).
Figure 7: Polarization curve.
Application Library path: Fuel_Cell_and_Electrolyzer_Module/Electrolyzers/soec_thermodynamics
Modeling Instructions
This tutorial models the current distribution in a solid oxide electrolyzer. The tutorial comprises two major parts. First, the thermodynamics and the electrode reactions are defined to provide the equilibrium potentials of the anode and the cathode, and a secondary (not concentration dependent) current distribution is modeled. In the second part, mass and momentum transport are added to model a concentration-dependent current distribution of the cell, where the mixture properties of the anode gas depends on the mole fractions of water and hydrogen.
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 > Primary and Secondary Current Distribution > Secondary Current Distribution (cd).
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
The model geometry is available as a parameterized geometry sequence in a separate MPH file. If you want to build it from scratch, follow the instructions in the section Appendix —Geometry Modeling Instructions. Otherwise load it from file using the following steps.
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 soec_thermodynamics_geom_sequence.mph.
3
In the Geometry toolbar, click  Build All.
4
In the Model Builder window, under Component 1 (comp1) click Geometry 1.
Global Definitions
Geometry Parameters
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
Some parameters were imported with the geometry sequence. Import some additional physics parameters from a text file.
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 soec_thermodynamics_physics_parameters.txt.
Add a thermodynamic system for the hydrogen/water gas mixture.
5
In the Physics toolbar, click  Thermodynamics and choose Thermodynamic System.
Select System
1
Go to the Select System window.
2
Click the Next button in the window toolbar.
Select Species
1
Go to the Select Species window.
2
In the Species list box, select water (7732-18-5, H2O).
3
Click  Add Selected.
4
In the Species list box, select hydrogen (1333-74-0, H2).
5
Click  Add Selected.
6
Click the Next button in the window toolbar.
Select Thermodynamic Model
1
Go to the Select Thermodynamic Model window.
2
Click the Finish button in the window toolbar.
Global Definitions
Gas System - H2 and H2O
1
In the Settings window for Thermodynamic System, type Gas System - H2 and H2O in the Label text field.
2
Right-click Global Definitions > Thermodynamics > Gas System - H2 and H2O and choose Generate Chemistry.
Select Species
1
Go to the Select Species window.
2
In the list, choose hydrogen and water.
3
Click  Add Selected.
4
Click the Next button in the window toolbar.
Chemistry Settings
1
Go to the Chemistry Settings window.
2
From the Mass transfer list, choose Concentrated species.
3
Click the Finish button in the window toolbar.
Definitions
Add an average operator which will be used later while plotting the polarization plot.
Average 1 (aveop1)
1
In the Definitions toolbar, click  Nonlocal Couplings and choose Average.
2
In the Settings window for Average, type aveop_an in the Operator name text field.
3
Locate the Source Selection section. From the Geometric entity level list, choose Boundary.
4
From the Selection list, choose Anode Current Collector.
Chemistry - H2 and H2O
1
In the Model Builder window, under Component 1 (comp1) click Chemistry (chem).
2
In the Settings window for Chemistry, type Chemistry - H2 and H2O in the Label text field.
3
Locate the Model Input section. Select the Enable electrode reactions checkbox.
4
From the E list, choose Electrode potential (cd).
5
Click to expand the Calculate Transport Properties section.
Electrode Reaction 1
Add the electrode reaction for the cathode. Use the (ads) annotation in the reaction formula to define that the oxygen ions do not belong the main (gas) phase.
1
In the Physics toolbar, click  Domains and choose Electrode Reaction.
2
In the Settings window for Electrode Reaction, locate the Reaction Formula section.
3
In the Formula (written as reduction) text field, type H2O+2e<=>H2+O(ads).
4
Click Apply.
5
Locate the Equilibrium Potential section. From the list, choose Automatic.
6
Locate the Electrode Kinetics section. In the i0,ref(T) text field, type i0_H2.
7
In the αa text field, type 1.5.
For the first computation we will assume constant mass fractions (concentrations) for the species.
8
In the Model Builder window, under Component 1 (comp1) click Chemistry - H2 and H2O (chem).
9
In the Settings window for Chemistry, locate the Species Matching section.
10
Find the Bulk species subsection. In the table, enter the following settings:
Species
Type
Mass fraction
Value (1)
From Thermodynamics
H2
Variable, from Reaction
User defined
0.5
H2
H2O
Variable, from Reaction
User defined
0.5
H2O
11
Find the Surface species subsection. In the table, enter the following settings:
Species
Type
Surface concentration
Value (mol/m^2)
O(ads)
Free species
User defined
1e-5[mol/m^2]
Global Definitions
Now add a second thermodynamics system for the anode oxygen gas.
In the Physics toolbar, click  Thermodynamics and choose Thermodynamic System.
Select System
1
Go to the Select System window.
2
Click the Next button in the window toolbar.
Select Species
1
Go to the Select Species window.
2
In the Species list box, select oxygen (7782-44-7, O2).
3
Click  Add Selected.
4
Click the Next button in the window toolbar.
Select Thermodynamic Model
1
Go to the Select Thermodynamic Model window.
2
Click the Finish button in the window toolbar.
Global Definitions
Gas System - O2
1
In the Settings window for Thermodynamic System, type Gas System - O2 in the Label text field.
2
Right-click Global Definitions > Thermodynamics > Gas System - O2 and choose Generate Chemistry.
Select Species
1
Go to the Select Species window.
2
In the list box, select oxygen.
3
Click  Add Selected.
4
Click the Next button in the window toolbar.
Chemistry Settings
1
Go to the Chemistry Settings window.
2
Click the Finish button in the window toolbar.
Chemistry - O2
1
In the Model Builder window, under Component 1 (comp1) click Chemistry 2 (chem2).
2
In the Settings window for Chemistry, type Chemistry - O2 in the Label text field.
3
Locate the Model Input section. Select the Enable electrode reactions checkbox.
4
From the E list, choose Electrode potential (cd).
Electrode Reaction 1
1
In the Physics toolbar, click  Domains and choose Electrode Reaction.
2
In the Settings window for Electrode Reaction, locate the Reaction Formula section.
3
In the Formula (written as reduction) text field, type O2+4e<=>2O(ads).
4
Click Apply.
5
Locate the Equilibrium Potential section. From the list, choose Automatic.
6
Locate the Electrode Kinetics section. In the i0,ref(T) text field, type i0_O2.
7
In the Model Builder window, under Component 1 (comp1) click Chemistry - O2 (chem2).
8
In the Settings window for Chemistry, locate the Species Matching section.
9
Find the Bulk species subsection. In the table, enter the following settings:
Species
Type
Mass fraction
Value (1)
From Thermodynamics
O2
Variable
User defined
1
O2
10
Find the Surface species subsection. In the table, enter the following settings:
Species
Type
Surface concentration
Value (mol/m^2)
O(ads)
Free species
User defined
1e-5[mol/m^2]
Secondary Current Distribution (cd)
Now define the current distribution in the gas diffusion electrodes and the electrolyte.
1
In the Model Builder window, under Component 1 (comp1) click Secondary Current Distribution (cd).
2
In the Settings window for Secondary Current Distribution, locate the Domain Selection section.
3
Click  Clear Selection.
4
Select Domains 1–3 only.
Electrolyte 1
1
In the Model Builder window, under Component 1 (comp1) > Secondary Current Distribution (cd) click Electrolyte 1.
2
In the Settings window for Electrolyte, locate the Electrolyte section.
3
From the σl list, choose User defined. In the associated text field, type sigma_l.
Porous Electrode - H2 and H2O (Cathode)
1
In the Physics toolbar, click  Domains and choose Porous Electrode.
2
In the Settings window for Porous Electrode, type Porous Electrode - H2 and H2O (Cathode) in the Label text field.
3
Locate the Domain Selection section. From the Selection list, choose Cathode.
4
Locate the Electrolyte Current Conduction section. From the σl list, choose User defined. In the associated text field, type sigma_l.
5
In the εl text field, type por_l.
6
Locate the Electrode Current Conduction section. From the σs list, choose User defined. In the associated text field, type sigma_s.
7
In the εs text field, type 1-por_l.
Porous Electrode Reaction 1
Couple the electrode reaction current density to the chemistry interface for the cathode as follows.
1
In the Model Builder window, click Porous Electrode Reaction 1.
2
In the Settings window for Porous Electrode Reaction, locate the Equilibrium Potential section.
3
In the Eeq text field, type chem.Eeq_er1.
4
Locate the Electrode Kinetics section. From the iloc,expr list, choose User defined. In the associated text field, type chem.iloc_er1.
5
Locate the Active Specific Surface Area section. In the av text field, type S.
Porous Electrode - O2 (Anode)
1
In the Model Builder window, right-click Porous Electrode - H2 and H2O (Cathode) and choose Duplicate.
2
In the Settings window for Porous Electrode, type Porous Electrode - O2 (Anode) in the Label text field.
3
Locate the Domain Selection section. From the Selection list, choose Anode.
Porous Electrode Reaction 1
1
In the Model Builder window, expand the Porous Electrode - O2 (Anode) node, then click Porous Electrode Reaction 1.
2
In the Settings window for Porous Electrode Reaction, locate the Equilibrium Potential section.
3
In the Eeq text field, type chem2.Eeq_er1.
4
Locate the Electrode Kinetics section. In the iloc,expr text field, type chem2.iloc_er1.
Electric Ground 1
1
In the Physics toolbar, click  Boundaries and choose Electric Ground.
2
In the Settings window for Electric Ground, locate the Boundary Selection section.
3
From the Selection list, choose Cathode Current Collector.
4
Locate the Contact Resistance section. Select the Include contact resistance checkbox.
5
In the Rc text field, type Rc.
Electrode Current 1
1
In the Physics toolbar, click  Boundaries and choose Electrode Current.
2
In the Settings window for Electrode Current, locate the Boundary Selection section.
3
From the Selection list, choose Anode Current Collector.
4
Locate the Electrode Current section. From the list, choose Average current density.
5
In the is,average text field, type I_avg.
6
Locate the Contact Resistance section. Select the Include contact resistance checkbox.
7
In the Rc text field, type Rc.
Global Definitions
Default Model Inputs
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
The physics settings for the first part of the tutorial are now complete. Add a user-defined mesh.
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, click to expand the Element Size Parameters section.
3
Locate the Element Size section. Click the Custom button.
4
Locate the Element Size Parameters section.
5
Select the Maximum element size checkbox. In the associated text field, type H_ch*0.8.
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
Select Domains 5–7 only.
5
Click  Build Selected.
6
Click  Build Selected.
Free Tetrahedral 1
1
In the Mesh toolbar, click  Free Tetrahedral.
2
In the Settings window for Free Tetrahedral, locate the Domain Selection section.
3
From the Geometric entity level list, choose Domain.
4
Select Domains 4 and 8 only.
5
Click  Build Selected.
Boundary Layers 1
Also add a boundary layer mesh at this stage. These are actually not needed for the first calculation, but will improve the accuracy and convergence of the solution for the second part of the tutorial when mass transport and convection has been added.
1
In the Mesh toolbar, click  Boundary Layers.
2
In the Settings window for Boundary Layers, locate the Geometric Entity Selection section.
3
From the Geometric entity level list, choose Domain.
4
From the Selection list, choose Channel Domains.
Boundary Layer Properties
1
In the Model Builder window, click Boundary Layer Properties.
2
In the Settings window for Boundary Layer Properties, locate the Boundary Selection section.
3
From the Selection list, choose Boundary Layer Boundaries.
4
Locate the Layers section. In the Number of layers text field, type 2.
5
In the Stretching factor text field, type 1.3.
6
From the Thickness specification list, choose First layer.
7
In the Thickness text field, type H_ch/10.
8
Click  Build Selected.
Free Triangular 1
1
In the Mesh toolbar, click  More Generators and choose Free Triangular.
2
In the Settings window for Free Triangular, locate the Boundary Selection section.
3
From the Selection list, choose Cathode Current Collector.
4
Click  Build Selected.
Swept 2
In the Mesh toolbar, click  Swept.
Distribution 1
1
Right-click Swept 2 and choose Distribution.
2
In the Settings window for Distribution, locate the Domain Selection section.
3
Click  Clear Selection.
4
Select Domain 3 only.
5
Locate the Distribution section. From the Distribution type list, choose Predefined.
6
In the Element ratio text field, type 2.
7
Select the Reverse direction checkbox.
Distribution 2
1
In the Model Builder window, right-click Swept 2 and choose Distribution.
2
In the Settings window for Distribution, locate the Domain Selection section.
3
Click  Clear Selection.
4
Select Domain 2 only.
5
Locate the Distribution section. In the Number of elements text field, type 2.
Distribution 3
1
In the Model Builder window, under Component 1 (comp1) > Mesh 1 > Swept 2 right-click Distribution 1 and choose Duplicate.
2
In the Settings window for Distribution, locate the Domain Selection section.
3
Click  Clear Selection.
4
Select Domain 1 only.
5
Locate the Distribution section. Clear the Reverse direction checkbox.
6
Click  Build All.
Study 1
The problem is now ready for solving.
1
In the Study toolbar, click  Compute.
Results
Electrode Potential with Respect to Ground (cd)
1
In the Model Builder window, under Results click Electrode Potential with Respect to Ground (cd).
2
In the Electrode Potential with Respect to Ground (cd) toolbar, click  Plot.
Inspect the potential plot. The plot should look as follows:
Component 1 (comp1)
Now start with the second part of the tutorial which adds a reacting flow interface to the cathode side of the model, and couples the distribution of hydrogen and water vapor to the electrochemistry.
Add Physics
1
In the Home toolbar, click  Add Physics to open the Add Physics window.
2
Go to the Add Physics window.
3
In the tree, select Chemical Species Transport > Reacting Flow in Porous Media > Transport of Concentrated Species.
4
Click the Add to Component 1 button in the window toolbar.
5
In the Home toolbar, click  Add Physics to close the Add Physics window.
Multiphysics
Reacting Flow 1 (nirf1)
1
In the Settings window for Reacting Flow, locate the Temperature section.
2
In the T text field, type T.
Materials
Note that a Porous Material node has been created automatically when adding an entry from the Reacting Flow in Porous Media branch of the Add Physics window.
Porous Material 1 (pmat1)
1
In the Model Builder window, expand the Component 1 (comp1) > Materials node, then click Porous Material 1 (pmat1).
2
In the Settings window for Porous Material, locate the Geometric Entity Selection section.
3
From the Selection list, choose Cathode.
4
Locate the Porosity section. In the εp text field, type por.
5
Locate the Homogenized Properties section. In the table, enter the following settings:
Property
Variable
Value
Unit
Property group
Permeability
kappa_iso ; kappaii = kappa_iso, kappaij = 0
kappa
m²
Basic
Brinkman Equations (br)
1
In the Model Builder window, under Component 1 (comp1) click Brinkman Equations (br).
2
Select Domains 3–8 only.
3
In the Settings window for Brinkman Equations, locate the Domain Selection section.
4
Click  Create Selection.
5
In the Create Selection dialog, type Gas domains in the Selection name text field.
6
Click OK.
7
In the Settings window for Brinkman Equations, locate the Physical Model section.
8
Clear the Neglect inertial term (Stokes flow) checkbox.
Fluid Properties 1
1
In the Physics toolbar, click  Domains and choose Fluid Properties.
2
In the Settings window for Fluid Properties, locate the Domain Selection section.
3
From the Selection list, choose Channel Domains.
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 Inlet.
4
Locate the Boundary Condition section. From the list, choose Mass flow.
5
Locate the Mass Flow section. In the m text field, type Mflux_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 Outlet.
4
Locate the Pressure Conditions section. Select the Normal flow checkbox.
Transport of Concentrated Species in Porous Media (tcs)
1
In the Model Builder window, under Component 1 (comp1) click Transport of Concentrated Species in Porous Media (tcs).
2
In the Settings window for Transport of Concentrated Species in Porous Media, locate the Domain Selection section.
3
From the Selection list, choose Gas domains.
4
Locate the Transport Mechanisms section. From the Diffusion model list, choose Maxwell–Stefan.
5
Click to expand the Dependent Variables section. In the Mass fractions (1) table, enter the following settings:
wH2
wH2O
Species Molar Masses 1
1
In the Model Builder window, under Component 1 (comp1) > Transport of Concentrated Species in Porous Media (tcs) click Species Molar Masses 1.
2
In the Settings window for Species Molar Masses, locate the Molar Mass section.
3
From the MwH2 list, choose Molar mass (chem/H2).
4
From the MwH2O list, choose Molar mass (chem/H2O).
Fluid 1
1
In the Model Builder window, under Component 1 (comp1) > Transport of Concentrated Species in Porous Media (tcs) > Porous Medium 1 click Fluid 1.
2
In the Settings window for Fluid, locate the Diffusion section.
3
In the table, enter the following settings:
Species 1
Species 2
Diffusivity
Diffusion coefficient (m^2/s)
wH2
wH2O
Maxwell-Stefan diffusivity, H2-H2O (chem)
comp1.chem.D_H2_H2O
4
From the Effective diffusivity model list, choose Bruggeman model.
5
Locate the Pore-Wall Interaction section. Select the Include pore-wall interaction checkbox.
6
In the dpore text field, type d_pore.
Porous Electrode Coupling 1
1
In the Physics toolbar, click  Domains and choose Porous Electrode Coupling.
2
In the Settings window for Porous Electrode Coupling, locate the Domain Selection section.
3
From the Selection list, choose Cathode.
Reaction Coefficients 1
1
In the Model Builder window, click Reaction Coefficients 1.
2
In the Settings window for Reaction Coefficients, locate the Reaction Current Source section.
3
From the iv list, choose Local current source, Porous Electrode Reaction 1 (cd/pce1/per1).
4
Locate the Stoichiometric Coefficients section. In the n text field, type 2.
5
In the νwH2 text field, type 1.
6
In the νwH2O text field, type -1.
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 Inlet.
4
Locate the Inflow section. From the Mixture specification list, choose Mass flow rates.
5
In the Jin,wH2O text field, type Mflux_in.
6
Click the  Show More Options button in the Model Builder toolbar.
7
In the Show More Options dialog, select Physics > Advanced Physics Options in the tree.
8
In the tree, select the checkbox for the node Physics > Advanced Physics Options.
9
Click OK.
10
In the Settings window for Inflow, click to expand the Boundary Mass Fraction Initial Values section.
11
In the ω0,bnd,wH2O text field, type 0.99.
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 Outlet.
Fluid 1
1
In the Physics toolbar, click  Domains and choose Fluid.
2
In the Settings window for Fluid, locate the Domain Selection section.
3
From the Selection list, choose Channel Domains.
4
Locate the Diffusion section. In the table, enter the following settings:
Species 1
Species 2
Diffusivity
Diffusion coefficient (m^2/s)
wH2
wH2O
Maxwell-Stefan diffusivity, H2-H2O (chem)
comp1.chem.D_H2_H2O
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
From the Mixture specification list, choose Mole fractions.
4
In the x0,wH2O text field, type 0.99.
Chemistry - H2 and H2O (chem)
1
In the Model Builder window, under Component 1 (comp1) click Chemistry - H2 and H2O (chem).
2
In the Settings window for Chemistry, locate the Species Matching section.
3
Find the Bulk species subsection. From the Species solved for list, choose Transport of Concentrated Species in Porous Media.
4
In the table, enter the following settings:
Species
Type
Mass fraction
Value (1)
From Thermodynamics
H2
Variable
wH2
Solved for
H2
H2O
Variable
wH2O
Solved for
H2O
Study 1
The concentration-dependent model is now ready for solving. Use a sequence of study steps, solving for the secondary current distribution first, then the flow, and finally the fully coupled problem. By solving for only one set of physics at a time in individual steps, suitable initial values automatically propagate to the final study step where the complete problem is solved.
Step 2: Stationary
1
In the Model Builder window, under Study 1 click Step 2: 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 checkboxes for Brinkman Equations (br) and Transport of Concentrated Species in Porous Media (tcs).
4
In the Solve for column of the table, under Component 1 (comp1) > Multiphysics, clear the checkbox for Reacting Flow 1 (nirf1).
Step 3: Stationary 2
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 checkboxes for Secondary Current Distribution (cd) and Transport of Concentrated Species in Porous Media (tcs).
4
In the Solve for column of the table, under Component 1 (comp1) > Multiphysics, clear the checkbox for Reacting Flow 1 (nirf1).
Step 4: Stationary 3
1
In the Study toolbar, click  Stationary.
Use Auxiliary sweep to generate the polarization plot.
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
I_avg (Average current density)
I_avg_max*range(0.1,0.1,1)
A/m^2
Solver Configurations
Remove the old study sequence and generate a new one.
In the Model Builder window, under Study 1 right-click Solver Configurations and choose Delete Configurations.
Solution 1 (sol1)
1
In the Study toolbar, click  Show Default Solver.
2
In the Model Builder window, expand the Solution 1 (sol1) node.
3
In the Model Builder window, expand the Study 1 > Solver Configurations > Solution 1 (sol1) > Stationary Solver 4 node.
4
Right-click Stationary Solver 4 and choose Fully Coupled.
5
In the Study toolbar, click  Compute.
The problem should solve in about two minutes.
Results
Start the postprocessing of the solution by inspecting and polishing the default plot for the velocity field.
Multislice 1
1
In the Model Builder window, expand the Results > Velocity (br) node, then click Multislice 1.
2
In the Settings window for Multislice, locate the Multiplane Data section.
3
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 H_cell-H_ch/2.
7
In the Velocity (br) toolbar, click  Plot.
Streamline 1
1
In the Model Builder window, right-click Velocity (br) and choose Streamline.
2
In the Settings window for Streamline, click Replace Expression in the upper-right corner of the Expression section. From the menu, choose Component 1 (comp1) > Brinkman Equations > Velocity and pressure > u,v,w - Velocity field.
3
Locate the Selection section. From the Selection list, choose Inlet.
4
Locate the Coloring and Style section. Find the Point style subsection. From the Type list, choose Arrow.
5
From the Arrow distribution list, choose Equal inverse time.
6
From the Color list, choose Black.
7
In the Velocity (br) toolbar, click  Plot.
Density
The density of the gas mixture will change as water is replaced by hydrogen in the gas stream. Plot the density as follows:
1
In the Results toolbar, click  3D Plot Group.
2
In the Settings window for 3D Plot Group, type Density in the Label text field.
Surface 1
1
Right-click Density 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) > Brinkman Equations > Material properties > br.rho - Density - kg/m³.
3
In the Density toolbar, click  Plot.
Viscosity
Also, the viscosity will change in the gas stream.
1
In the Results toolbar, click  3D Plot Group.
2
In the Settings window for 3D Plot Group, type Viscosity in the Label text field.
Surface 1
1
Right-click Viscosity 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) > Brinkman Equations > Material properties > br.mu - Dynamic viscosity - Pa·s.
3
In the Viscosity toolbar, click  Plot.
Concentration, H2, Surface (tcs)
Default plots are created for the hydrogen and water mole fractions.
Surface 1
1
In the Model Builder window, expand the Concentration, H2, Surface (tcs) 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 ConopiformisZero.
Surface 1
1
In the Model Builder window, expand the Concentration, H2O, Surface (tcs) 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 Prionace.
Mole Fraction and Flux, H2
Create a plot for the hydrogen mole fraction and flux as follows:
1
In the Results toolbar, click  3D Plot Group.
2
In the Settings window for 3D Plot Group, type Mole Fraction and Flux, H2 in the Label text field.
3
Locate the Plot Settings section. Clear the Plot dataset edges checkbox.
Streamline 1
1
Right-click Mole Fraction and Flux, H2 and choose Streamline.
2
In the Settings window for Streamline, click Replace Expression in the upper-right corner of the Expression section. From the menu, choose Component 1 (comp1) > Transport of Concentrated Species in Porous Media > Species wH2 > Fluxes > tcs.tflux_wH2x,...,tcs.tflux_wH2z - Total flux.
3
Locate the Streamline Positioning section. In the Number text field, type 30.
4
Locate the Selection section. From the Selection list, choose Outlet.
5
Locate the Coloring and Style section. Find the Point style subsection. From the Type list, choose Arrow.
6
From the Arrow distribution list, choose Equal inverse time.
7
From the Color list, choose Gray.
Selection 1
1
Right-click Streamline 1 and choose Selection.
2
In the Settings window for Selection, locate the Selection section.
3
From the Selection list, choose Channel Domains.
Volume 1
1
In the Model Builder window, right-click Mole Fraction and Flux, H2 and choose Volume.
2
In the Settings window for Volume, click Replace Expression in the upper-right corner of the Expression section. From the menu, choose Component 1 (comp1) > Transport of Concentrated Species in Porous Media > Species wH2 > tcs.x_wH2 - Mole fraction - 1.
3
Locate the Coloring and Style section. From the Color table list, choose ConopiformisZero.
Selection 1
1
Right-click Volume 1 and choose Selection.
2
In the Settings window for Selection, locate the Selection section.
3
From the Selection list, choose Cathode.
4
Click the  Show Grid button in the Graphics toolbar.
5
In the Mole Fraction and Flux, H2 toolbar, click  Plot.
Cross-Sectional Electrolyte Current Density
Plot the current distribution across the electrolyte layer as follows:
1
In the Results toolbar, click  3D Plot Group.
2
In the Settings window for 3D Plot Group, type Cross-Sectional Electrolyte Current Density in the Label text field.
3
Locate the Plot Settings section. Clear the Plot dataset edges checkbox.
Slice 1
1
Right-click Cross-Sectional Electrolyte Current Density and choose Slice.
2
In the Settings window for Slice, click Replace Expression in the upper-right corner of the Expression section. From the menu, choose Component 1 (comp1) > Secondary Current Distribution > Electrolyte current density vector - A/m² > cd.Ilz - Electrolyte current density vector, z-component.
3
Locate the Plane Data section. From the Plane list, choose xy-planes.
4
From the Entry method list, choose Coordinates.
5
In the z-coordinates text field, type H_gde+H_el/2.
6
Locate the Coloring and Style section. From the Color table list, choose MetasepiaBlue.
7
In the Cross-Sectional Electrolyte Current Density toolbar, click  Plot.
Polarization Plot
Finally, add the Polarization Plot.
1
In the Results toolbar, click  1D Plot Group.
2
In the Settings window for 1D Plot Group, type Polarization Plot in the Label text field.
3
Click to expand the Title section. From the Title type list, choose Manual.
4
In the Title text area, type Polarization plot.
5
Locate the Plot Settings section. Select the x-axis label checkbox.
6
Select the y-axis label checkbox.
7
In the x-axis label text field, type Average current density (A/cm<sup>2</sup>).
8
In the y-axis label text field, type Cell voltage (V).
Global 1
1
In the Polarization Plot toolbar, click  Global.
2
In the Settings window for Global, click Add Expression in the upper-right corner of the y-Axis Data section. From the menu, choose Component 1 (comp1) > Secondary Current Distribution > cd.phis0_ec1 - Electric potential on boundary - V.
3
Locate the x-Axis Data section. From the Parameter list, choose Expression.
4
In the Expression text field, type -aveop_an(cd.nIs).
5
In the Unit field, type A/cm^2.
6
Click to expand the Legends section. Clear the Show legends checkbox.
7
In the Polarization Plot toolbar, click  Plot.
Follow the commands below to improve some plot appearances and remove redundant plots.
Streamline 1
1
In the Model Builder window, expand the Concentration, H2, Streamline (tcs) node, then click Streamline 1.
2
In the Settings window for Streamline, locate the Streamline Positioning section.
3
In the Points text field, type 200.
4
Locate the Coloring and Style section. Find the Line style subsection. From the Type list, choose Tube.
5
Select the Radius scale factor checkbox. In the associated text field, type 2e-5.
6
Find the Point style subsection. From the Type list, choose None.
Color Expression
1
In the Model Builder window, expand the Streamline 1 node, then click Color Expression.
2
In the Settings window for Color Expression, locate the Coloring and Style section.
3
From the Color table list, choose ConopiformisZero.
4
In the Concentration, H2, Streamline (tcs) toolbar, click  Plot.
Streamline 1
1
In the Model Builder window, expand the Concentration, H2O, Streamline (tcs) node, then click Streamline 1.
2
In the Settings window for Streamline, locate the Streamline Positioning section.
3
In the Points text field, type 200.
4
Locate the Coloring and Style section. Find the Line style subsection. From the Type list, choose Tube.
5
Select the Radius scale factor checkbox. In the associated text field, type 2e-5.
6
Find the Point style subsection. From the Type list, choose None.
Color Expression
1
In the Model Builder window, expand the Streamline 1 node, then click Color Expression.
2
In the Settings window for Color Expression, locate the Coloring and Style section.
3
From the Color table list, choose Prionace.
4
In the Concentration, H2O, Streamline (tcs) toolbar, click  Plot.
Electrode Current Density (cd), Electrolyte Current Density (cd), Electrolyte Potential (cd)
1
In the Model Builder window, under Results, Ctrl-click to select Electrolyte Potential (cd), Electrolyte Current Density (cd), and Electrode Current Density (cd).
2
Right-click and choose Delete.
Appendix —Geometry 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
Click  Done.
Global Definitions
Geometry Parameters
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.
3
Locate the Parameters section. Click  Load from File.
4
Browse to the model’s Application Libraries folder and double-click the file soec_thermodynamics_geom_parameters.txt.
Geometry 1
Anode
1
In the Geometry toolbar, click  Block.
2
In the Settings window for Block, type Anode in the Label text field.
3
Locate the Size and Shape section. In the Width text field, type W_cell.
4
In the Depth text field, type D_cell.
5
In the Height text field, type H_gde.
6
Locate the Selections of Resulting Entities section. Select the Resulting objects selection checkbox.
Electrolyte
1
Right-click Anode and choose Duplicate.
2
In the Settings window for Block, type Electrolyte in the Label text field.
3
Locate the Position section. In the z text field, type H_gde.
Cathode
1
Right-click Electrolyte and choose Duplicate.
2
In the Settings window for Block, type Cathode in the Label text field.
3
Locate the Position section. In the z text field, type H_gde+H_el.
Work Plane 1 (wp1)
1
In the Geometry toolbar, click  Work Plane.
2
In the Settings window for Work Plane, locate the Plane Definition section.
3
In the z-coordinate text field, type H_cell-H_ch.
Work Plane 1 (wp1) > Plane Geometry
In the Model Builder window, click Plane Geometry.
Work Plane 1 (wp1) > Rectangle 1 (r1)
1
In the Work Plane toolbar, click  Rectangle.
2
In the Settings window for Rectangle, locate the Size and Shape section.
3
In the Width text field, type W_ch.
4
In the Height text field, type N_ch*(W_ch+W_rib).
5
Locate the Position section. In the xw text field, type W_rib/2.
6
In the yw text field, type W_rib/2.
Work Plane 1 (wp1) > Rectangle 2 (r2)
1
Right-click Component 1 (comp1) > Geometry 1 > Work Plane 1 (wp1) > Plane Geometry > Rectangle 1 (r1) and choose Duplicate.
2
In the Settings window for Rectangle, locate the Position section.
3
In the xw text field, type W_rib/2+L_ch+W_ch.
4
In the yw text field, type -W_rib/2.
Work Plane 1 (wp1) > Rectangle 3 (r3)
1
In the Work Plane toolbar, click  Rectangle.
2
In the Settings window for Rectangle, locate the Size and Shape section.
3
In the Width text field, type L_ch.
4
In the Height text field, type W_ch.
5
Locate the Position section. In the xw text field, type W_rib/2+W_ch.
6
In the yw text field, type W_rib/2.
Work Plane 1 (wp1) > Array 1 (arr1)
1
In the Work Plane toolbar, click  Transforms and choose Array.
2
Select the object r3 only.
3
In the Settings window for Array, locate the Size section.
4
In the yw size text field, type N_ch.
5
Locate the Displacement section. In the yw text field, type W_ch+W_rib.
6
In the Work Plane toolbar, click  Build All.
Channel Domains
1
In the Model Builder window, right-click Geometry 1 and choose Extrude.
2
In the Settings window for Extrude, type Channel Domains in the Label text field.
3
Locate the Distances section. In the table, enter the following settings:
Distances (m)
H_ch
4
Locate the Selections of Resulting Entities section. Select the Resulting objects selection checkbox.
Form Union (fin)
In the Model Builder window, right-click Form Union (fin) and choose Build Selected.
Inlet
1
In the Geometry toolbar, click  Selections and choose Explicit Selection.
2
In the Settings window for Explicit Selection, type Inlet in the Label text field.
3
Locate the Entities to Select section. From the Geometric entity level list, choose Boundary.
4
On the object fin, select Boundary 19 only.
Outlet
1
In the Geometry toolbar, click  Selections and choose Explicit Selection.
2
In the Settings window for Explicit Selection, type Outlet in the Label text field.
3
Locate the Entities to Select section. From the Geometric entity level list, choose Boundary.
4
On the object fin, select Boundary 42 only.
Cathode Current Collector
1
In the Geometry toolbar, click  Selections and choose Explicit Selection.
2
In the Settings window for Explicit Selection, type Cathode Current Collector in the Label text field.
3
Locate the Entities to Select section. From the Geometric entity level list, choose Boundary.
4
On the object fin, select Boundaries 10, 26, 33, and 40 only.
Anode Current Collector
1
In the Geometry toolbar, click  Selections and choose Explicit Selection.
2
In the Settings window for Explicit Selection, type Anode Current Collector in the Label text field.
3
Locate the Entities to Select section. From the Geometric entity level list, choose Boundary.
4
On the object fin, select Boundary 3 only.
Channel Domain Boundaries
1
In the Geometry toolbar, click  Selections and choose Adjacent Selection.
2
In the Settings window for Adjacent Selection, locate the Input Entities section.
3
Click  Add.
4
In the Add dialog, select Channel Domains in the Input selections list.
5
Click OK.
6
In the Settings window for Adjacent Selection, type Channel Domain Boundaries in the Label text field.
Boundary Layer Boundaries
1
In the Geometry toolbar, click  Selections and choose Difference Selection.
2
In the Settings window for Difference Selection, type Boundary Layer Boundaries in the Label text field.
3
Locate the Geometric Entity Level section. From the Level list, choose Boundary.
4
Locate the Input Entities section. Click the  Add button for Selections to add.
5
In the Add dialog, select Channel Domain Boundaries in the Selections to add list.
6
Click OK.
7
In the Settings window for Difference Selection, locate the Input Entities section.
8
Click the  Add button for Selections to subtract.
9
In the Add dialog, in the Selections to subtract list, choose Inlet and Outlet.
10
Click OK.