Models - Cfd.pem Electrolyzer

Created in COMSOL Multiphysics 5.
Polymer Electrolyte Membrane Electrolyzer
This model is licensed under the COMSOL Software License Agreement 5.6.
All trademarks are the property of their respective owners. See www.comsol.com/trademarks.
Introduction
In a polymer electrolyte membrane electrolyzer cell (PEMEC), the two electrode
compartments are separated by a polymer membrane, coated by porous gas diffusion
electrodes. Liquid water is fed to the anode side, forming oxygen gas on the anode side
and hydrogen gas on the cathode side.
The respective designs of the flow field patterns are important in order to obtain a uniform
distribution of flow, in combination with low pressure drops, during operation.
In this example, the mixture model is used to model the two-phase fluid dynamics on the
anode side of a PEMEC.
The model geometry and operating condition were taken from Ref. 1, with added gravity
and zero tangential (no slip) conditions for all channel walls. The single-phase results for
a 60 ml/min flow rate were verified versus Ref. 2.
Model Definition
Figure 1 shows the model geometry. From the circular inlet (located at the top boundary
of the cylinder), liquid water is led into a manifold, which in turn distributes the flow over
23 channels. Oxygen gas is produced at the anode electrode, located below the
2 | POLYMER ELECTROLYTE MEMBRANE ELECTROLYZER

23 channels. The two-phase oxygen gas/liquid water mixture exits the cell through the
exit manifold.
Outlet (liquid water + oxygen gas)
Inlet (liquid water)
Figure 1: Geometry of the anode flow field.
The model is set up using the Mixture Model, Laminar Flow interface, with liquid water
defining the continuous phase and the oxygen gas bubbles the dispersed phase.
Incompressible and isothermal conditions are assumed.
The inlet flow rate of liquid water is 260 ml/min. This is defined using an Inlet boundary
condition.
At the electrode surface/channel boundaries liquid water is consumed, and oxygen gas is
produced according to
+ -
2H 2 O(l) → O 2 (g) + 4H + 4e
The protons are transported through the polymer membrane, dividing the two electrode
compartments, over to the cathode side of the electrolyzer cell. In addition to the oxygen
gas produced, there is hence a net mass outflux due to the proton transport at the
electrode surface/channel boundaries.
The combined oxygen gas production/total mass outflux is defined using an Inlet
boundary condition node, assuming a total oxygen production of 5 mg/s.

A pressure condition is used at the outlet boundaries. No-slip wall conditions are set for
all other boundaries.
Finally, the buoyancy effects due to gravity are included in the model using a Gravity node,
with the gravity vector pointing downward in the z direction.
The model is solved in two steps. First the single-phase (pure liquid water, no oxygen
production) stationary flow is computed using a stationary solver. This solution is then
used as initial conditions for a 10 s time-dependent simulation, where the oxygen
production is ramped up to full production from 0 during the first second.
Results and Discussion

Figure 2 shows a slice plot of the mass-averaged velocity magnitude.
Figure 2: Slice plot of the velocity magnitude at t = 10 s.
The highest velocities are found in the inlet/outlet manifold channels. A good practice
when assuming laminar flow is to check the Reynolds number of the computed results.
The Reynolds number Re is defined as
ρuD
Re= ------------
μ

where ρ is the density, u the velocity, μ the dynamic viscosity and D the characteristic
length.
Given that the width of the channels, 2 mm, is larger than the height, 0.889 mm, we
choose the doubled height (an approximation valid for a wide duct) as the characteristic
length D. At the inlet, where we have pure water and the density-to-dynamic viscosity is
the highest, the maximum velocity is about 1.3 m/s.
The Reynolds number for the inlet manifold channels becomes (all parameter values using
the corresponding SI units)
3 –3
ρuD 10 × 1.3 × ( 2 × 0.889 × 10 )-
Re= ------------ ≈ ----------------------------------------------------------------------------- ≈ 2300
μ 10
–3
A similar calculation for the outlet renders lower values. Reynolds numbers in the range of
2300 and lower, indicates that turbulence should not have to be considered for the given
geometry and flow rates.
Figure 3 shows the gas volume fraction due to evolved oxygen in the cell at 10 s.
Figure 3: Gas volume fraction in the cell at t = 10 s.
The gas volume fraction approaches 100% at the end of the electrode flow channels located
at the middle of the flow field.

Figure 4 and Figure 5 show the pressure drop in the anode flow field at t = 0 and t = 10 s,
respectively. The pressure drop is significantly reduced as a result of the oxygen gas
production.
Figure 4: Pressure drop in the cell at t = 0 s.

Figure 5: Pressure drop in the cell at t = 10 s.
Finally, Figure 6 shows the velocity magnitudes at half the length (y direction) and half in
the height (z direction) of the electrode channels at various times. This plot is important
since it indicates the uniformness of the flow distribution over the individual channels. As
can be seen, the flow distribution not particularly uniform for pure water (t=0 s), but gets
significantly even less uniform when the gas production starts (t=1 and 2 s). At t=10 s the
distribution has relaxed back to a somewhat more uniform profile, but still less uniform
than for pure water. It is also seen that the flow field distribution does not change

significantly between 2 s and 10 s. This indicates that a stationary flow distribution is
established fairly soon after full oxygen production has been reached at 1 s.
Figure 6: Individual channel velocities at various times.
Notes About the COMSOL Implementation

The local oxygen flux is multiplied by a smoothed step function, going from 1 to 0 when
the volume fraction of oxygen approaches 1. The smoothing improves convergence.
Reference
1. J. Nie and Y. Chen, “Numerical modeling of three-dimensional two-phase gas-liquid
flow in the field plate of a PEM electrolysis cell,” Int. J. Hydrogen Energy, vol. 35,
pp. 3183–3197, 2010.
2. J. Nie, Y. Chen, S. Cohen, B. Carter, and R. Boehm, “Numerical and experimental

study of three-dimensional fluid flow in the bipolar plate of a PEM electrolysis cell,” Int.
J. Thermal Sciences, vol. 48, pp. 1914–1922, 2009.
Application Library path: CFD_Module/Multiphase_Flow/pem_electrolyzer

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 Fluid Flow>Multiphase Flow>Mixture Model>
Mixture Model, Laminar Flow (mm).
3 Click Add.
4 Click Study.
5 In the Select Study tree, select General Studies>Time Dependent.
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 with the
following steps.
1 In the Geometry toolbar, click Insert Sequence.

2 Browse to the model’s Application Libraries folder and double-click the file
pem_electrolyzer_geom_sequence.mph.

3 In the Geometry toolbar, click Build All.
GLOBAL DEFINITIONS
Use the parameterization to reduce the number of channels and the channel lengths when
setting up the model. It is often a good practice to start a modeling project on a reduced
geometry size (or dimension). This saves time and computational resources while
troubleshooting.
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 In the table, enter the following settings:
Name Expression Value Description

N_ch 3 3 Number of electrode channels
L_ch 118*h_a/5 0.02098 m Electrode channel lenghts
GEOMETRY 1
1 In the Geometry toolbar, click Build All.

2 Click the Zoom Extents button in the Graphics toolbar.
GLOBAL DEFINITIONS
Load some more physics parameters and variables from text files. Note that parameters and
variables defining the oxygen and water flows are scaled with the geometric parameters.
Parameters 1
3 Click Load from File.
pem_electrolyzer_parameters.txt.
DEFINITIONS
Variables 1
1 In the Model Builder window, under Component 1 (comp1) right-click Definitions and
choose Variables.
2 In the Settings window for Variables, locate the Variables section.

pem_electrolyzer_variables.txt.
The variables make use of a step and a ramp function. These have not yet been defined,
hence some of the loaded variable expressions are marked in orange. Define the missing
functions as follows:
Step 1 (step1)
1 In the Home toolbar, click Functions and choose Local>Step.
2 In the Settings window for Step, locate the Parameters section.
3 In the Location text field, type 0.95.
4 In the From text field, type 1.
5 In the To text field, type 0.
6 Click to expand the Smoothing section. In the Size of transition zone text field, type 0.1.
7 Click Plot.
The step function is used to set the oxygen flux to zero locally when the gas volume
fraction approaches 1. Smoothing is important in order to avoid discrete jumps in the
flux. We will decrease the smoothing later when solving for the full model.
Ramp 1 (rm1)
1 In the Home toolbar, click Functions and choose Local>Ramp.
2 In the Settings window for Ramp, locate the Parameters section.

3 Select the Cutoff check box.
4 Click Plot.
The ramp function is used to ramp up the oxygen flux from zero when the time-
dependent solver starts. This shortens the computational time.
DEFINITIONS
In the Model Builder window, collapse the Component 1 (comp1)>Definitions node.
ADD MATERIAL
1 In the Home toolbar, click Add Material to open the Add Material window.
Add liquid water and oxygen gas from the Material Library. Note that the order is
important - Add water first.
2 Go to the Add Material window.
3 In the tree, select Built-in>Water, liquid.
4 Click Add to Component 1 (comp1).
5 In the tree, select Liquids and Gases>Gases>Oxygen.
6 Click Add to Component 1 (comp1).

7 In the Home toolbar, click Add Material to close the Add Material window.
The oxygen node under Materials should now have a small warning symbol in the
Model Builder Tree. This is because the selection of this node is zero. This is expected
at this point.
Since the water node was added first, it got assigned to all domains by default.
MIXTURE MODEL, LAMINAR FLOW (MM)

1 In the Settings window for Mixture Model, Laminar Flow, locate the Physical Model
section.
2 From the Dispersed phase list, choose Liquid droplets/bubbles.
3 From the Slip model list, choose Schiller-Naumann.
In this model, water is the continuous phase, and oxygen the dispersed phase.
Mixture Properties 1
1 In the Model Builder window, under Component 1 (comp1)>Mixture Model,
Laminar Flow (mm) click Mixture Properties 1.
2 In the Settings window for Mixture Properties, locate the Materials section.
3 From the Continuous phase list, choose Water, liquid (mat1).
4 From the Dispersed phase list, choose Oxygen (mat2).
5 Locate the Dispersed Phase Properties section. From the ρd list, choose User defined. In
the associated text field, type rhoO2.
6 In the dd text field, type D_bubbles.
7 Locate the Mixture Model section. From the Mixture viscosity model list, choose
Volume averaged.
Inlet - Liquid Water

1 In the Physics toolbar, click Boundaries and choose Inlet.
2 In the Settings window for Inlet, type Inlet - Liquid Water in the Label text field.
3 Locate the Boundary Selection section. From the Selection list, choose Inlet.
4 Locate the Velocity section. In the J0 text field, type Flow_rate/(pi*R_in^2).
Inlet - Electrode Surface Oxygen Evolution

1 In the Physics toolbar, click Boundaries and choose Inlet.
2 In the Settings window for Inlet, type Inlet - Electrode Surface Oxygen
Evolution in the Label text field.
3 Locate the Boundary Selection section. From the Selection list, choose Electrode Surface.

4 Locate the Velocity section. In the J0 text field, type mixture_flow.
5 Locate the Dispersed Phase Boundary Condition section. From the
Dispersed phase boundary condition list, choose Dispersed phase flux.
6 In the Nφd text field, type disp_flow.
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.
Gravity 1
1 In the Physics toolbar, click Domains and choose Gravity.
The cell is oriented so that the z direction points upward.
2 In the Settings window for Gravity, locate the Domain Selection section.
3 From the Selection list, choose All domains.
4 Locate the Gravity section. Specify the g vector as
0 x
0 y
-g_const z
MESH 1
Manual meshing is required for a geometry of this complexity. Use a swept mesh along the
electrode channels, and free tetraheadal meshing for the remaining domains.
1 In the Model Builder window, under Component 1 (comp1) click Mesh 1.

2 In the Settings window for Mesh, locate the Mesh Settings section.
3 From the Sequence type list, choose Physics-controlled mesh.
Swept 1
Right-click Component 1 (comp1)>Mesh 1 and choose Swept.
Size
1 In the Settings window for Size, locate the Element Size section.
2 From the Predefined list, choose Finer.
Swept 1
1 In the Model Builder window, click Swept 1.

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 Channels Above Electrode Surface.
5 Click to expand the Source Faces section. From the Selection list, choose
Inlets to Electrode Channels .
6 Click to expand the Destination Faces section. From the Selection list, choose
Outlets from Electrode Channels.
Distribution 1
1 Right-click Swept 1 and choose Distribution.
2 In the Settings window for Distribution, locate the Domain Selection section.
4 Locate the Distribution section. In the Number of elements text field, type floor(L_ch/
(0.5*w_ch)).
Size 1
1 Right-click Swept 1 and choose Size.
3 Click the Custom button.
4 Locate the Geometric Entity Selection section. From the Geometric entity level list,
choose Boundary.
5 From the Selection list, choose Inlets to Electrode Channels .
6 Locate the Element Size Parameters section. Select the Maximum element size check box.
7 In the associated text field, type h_a/4.

8 Click Build Selected.
Free Tetrahedral 1
1 In the Model Builder window, right-click Mesh 1 and choose Free Tetrahedral.
2 In the Settings window for Free Tetrahedral, locate the Domain Selection section.
4 From the Selection list, choose Manifolds.
5 Click to expand the Scale Geometry section. In the z-direction scale text field, type 2.
Size 1
1 Right-click Free Tetrahedral 1 and choose Size.
5 In the associated text field, type w_ch/4.

Free Tetrahedral 2
1 In the Model Builder window, right-click Mesh 1 and choose Free Tetrahedral.
2 In the Settings window for Free Tetrahedral, locate the Scale Geometry section.
3 In the z-direction scale text field, type 2.
Size 1
1 Right-click Free Tetrahedral 2 and choose Size.
5 In the associated text field, type h_a.

Boundary Layers 1
Add boundary layer meshing to resolve steep velocity gradients at the inlet and outlet
regions to the electrode channels and along the walls..
1 In the Model Builder window, right-click Mesh 1 and choose Boundary Layers.
2 In the Settings window for Boundary Layers, locate the Geometric Entity Selection section.
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 Inlets and Outlets to Electrode Channels.
4 Locate the Boundary Layer Properties section. In the Number of boundary layers text
field, type 4.
5 From the Thickness of first layer list, choose Manual.
6 In the Thickness text field, type w_ch/15.

Boundary Layers 2
1 In the Model Builder window, right-click Mesh 1 and choose Boundary Layers.
2 In the Settings window for Boundary Layers, locate the Geometric Entity Selection section.
4 From the Selection list, choose Electrode Channels and Manifolds.
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 Layer Properties section.
3 In the Number of boundary layers text field, type 2.
4 From the Thickness of first layer list, choose Manual.
5 In the Thickness text field, type h_a/8.

6 Click Build All.
STUDY 1
The problem is now ready for solving. Add a Stationary study step to first solve for the
velocity and pressure fields for pure liquid water. This solution will then be used as initial
values for the time-dependent simulation.
Stationary
1 In the Study toolbar, click Study Steps and choose Stationary>Stationary.
2 Right-click Study 1>Step 2: Stationary and choose Move Up.
Step 2: Time Dependent

1 In the Model Builder window, click Step 2: Time Dependent.
2 In the Settings window for Time Dependent, locate the Study Settings section.
3 In the Output times text field, type 0 1 2 10.
Solution 1 (sol1)
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 1.
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 1 node, then click Volume fraction,
dispersed phase (comp1.phid).
6 In the Settings window for Field, locate the General section.
7 Clear the Solve for this field check box.
8 In the Model Builder window, click Squared slip velocity (comp1.slipvel).
By using manual scales for the volume fraction of the dispersed phase and the slip
velocity, the computation time can be reduced. (Suitable values for the slip velocity was
derived by solving the problem using automatic scales first, and then evaluating the
solution. This is another reason for why it makes sense to start with a small geometry
first.)
12 In the Settings window for Field, locate the Scaling section.
13 From the Method list, choose Manual.
17 In the Scale text field, type 1e-4.
18 In the Model Builder window, click Study 1.
19 In the Settings window for Study, locate the Study Settings section.
20 Clear the Generate default plots check box.
21 In the Study toolbar, click Compute.
The model should take about half an hour to solve.
RESULTS
Create plots for the pressure, velocity and gas volume fraction as follows:
Pressure
1 In the Home toolbar, click Add Plot Group and choose 3D Plot Group.

2 In the Settings window for 3D Plot Group, type Pressure in the Label text field.
Surface 1
1 Right-click Pressure 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)>Mixture Model,
Laminar Flow>VelAndPressure>p - Pressure - Pa.
3 In the Pressure toolbar, click Plot.
Velocity
2 In the Settings window for 3D Plot Group, type Velocity in the Label text field.
Slice 1
1 Right-click Velocity 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)>Mixture Model,
Laminar Flow>VelAndPressure>mm.U - Mass-averaged velocity field - m/s.
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_a/2.

6 In the Velocity toolbar, click Plot.
Gas Volume Fraction

2 In the Settings window for 3D Plot Group, type Gas Volume Fraction in the Label text
field.
Surface 1
1 Right-click Gas Volume Fraction and choose Surface.

2 In the Gas Volume Fraction toolbar, click Plot.
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_ch/2.
4 In row Point 1, set y to L_ch/2.
5 In row Point 1, set z to h_a/2.
6 In row Point 2, set x to N_ch*w_ch*2-3*w_ch/2.
7 In row Point 2, set y to L_ch/2.
8 In row Point 2, set z to h_a/2.
Velocity in Electrode Channels

1 In the Results toolbar, click 1D Plot Group.
2 In the Settings window for 1D Plot Group, type Velocity in Electrode Channels in
the Label text field.
3 Locate the Data section. From the Dataset list, choose Cut Line 3D 1.
Line Graph 1
1 Right-click Velocity in Electrode Channels and choose Line Graph.

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)>Mixture Model,
Laminar Flow>VelAndPressure>mm.U - Mass-averaged velocity field - m/s.
3 Locate the x-Axis Data section. From the Parameter list, choose Expression.
4 In the Expression text field, type x.
5 Click to expand the Legends section. Select the Show legends check box.

1 In the Model Builder window, click Velocity in Electrode Channels.
2 In the Settings window for 1D Plot Group, locate the Legend section.
3 From the Position list, choose Lower right.
4 In the Velocity in Electrode Channels toolbar, click Plot.
GLOBAL DEFINITIONS
Now model the full geometry. Go back and set the number of channels and the cell length
to their original values.
Parameters 1

Name Expression Value Description

N_ch 23 23 Number of electrode channels
L_ch 118*h_a 0.1049 m Electrode channel lenghts
GEOMETRY 1
1 In the Home toolbar, click Build All.
MESH 1
Inspect the mesh after the geometry change.

In the Model Builder window, under Component 1 (comp1) right-click Mesh 1 and choose
Build All.
DEFINITIONS
Step 1 (step1)
Decrease the smoothing of the step funciton.
1 In the Model Builder window, under Component 1 (comp1)>Definitions click

Step 1 (step1).
2 In the Settings window for Step, locate the Parameters section.
3 In the Location text field, type 0.995.
4 Locate the Smoothing section. In the Size of transition zone text field, type 0.01.

5 Click Plot.
STUDY 1
Solution 1 (sol1)
Before solving, make a copy of the solution for the small geometry to keep it for future
reference.
1 In the Model Builder window, under Study 1>Solver Configurations right-click

Solution 1 (sol1) and choose Solution>Copy.
Solution - Small Geometry

1 In the Model Builder window, under Study 1>Solver Configurations click Solution 1 -
Copy 1 (sol3).
2 In the Settings window for Solution, type Solution - Small Geometry in the Label
text field.
Solver Configurations
Reset the solver sequence in order to obtain the default solver for the new problem size.
Solution 1 (sol1)
1 In the Model Builder window, right-click Solver Configurations and choose
Reset Solver to Default.
2 In the Settings window for Dependent Variables, locate the General section.

3 From the Defined by study step list, choose User defined.
16 In the Scale text field, type 1e-5.
Step 2: Time Dependent

The problem will take several hours to solve. Plot the gas volume fraction while solving in
order to monitor the solution process.
1 In the Model Builder window, click Step 2: Time Dependent.

2 In the Settings window for Time Dependent, click to expand the Results While Solving
section.
3 Select the Plot check box.
4 From the Plot group list, choose Gas Volume Fraction.
5 From the Update at list, choose Time steps taken by solver.
The full problem is now ready for solving.
6 In the Home toolbar, click Compute.

RESULTS
Pressure
2 In the Settings window for 3D Plot Group, locate the Data section.
3 From the Time (s) list, choose 0.

5 From the Time (s) list, choose 10.
Velocity
1 In the Model Builder window, click Velocity.
2 In the Velocity toolbar, click Plot.

Gas Volume Fraction
1 In the Model Builder window, click Gas Volume Fraction.
2 In the Gas Volume Fraction toolbar, click Plot.

1 In the Model Builder window, click Velocity in Electrode Channels.

2 In the Velocity in Electrode Channels toolbar, click Plot.
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
Parameters 1
pem_electrolyzer_geom_sequence_parameters.txt.

GEOMETRY 1
Work Plane 1 (wp1)

In the Geometry toolbar, click Work Plane.
Work Plane 1 (wp1)>Plane Geometry

In the Model Builder window, click Plane Geometry.
Work Plane 1 (wp1)>Circle 1 (c1)

1 In the Work Plane toolbar, click Circle.
2 In the Settings window for Circle, locate the Size and Shape section.
3 In the Radius text field, type R_in.
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 L_inout*3/4.
4 In the Height text field, type w_ch.
5 Locate the Position section. In the yw text field, type -w_ch*2.5.

3 In the Width text field, type L_inout*1/4.
5 Locate the Position section. In the xw text field, type L_inout*3/4.
6 In the yw text field, type -w_ch*2.5.
Work Plane 1 (wp1)>Array 1 (arr1)

1 In the Work Plane toolbar, click Transforms and choose Array.
2 Select the objects r1 and r2 only.
3 In the Settings window for Array, locate the Size section.

4 In the yw size text field, type 3.
5 Locate the Displacement section. In the yw text field, type 2*w_ch.
Work Plane 1 (wp1)>Polygon 1 (pol1)

1 In the Work Plane toolbar, click Polygon.
2 In the Settings window for Polygon, locate the Coordinates section.
xw (m) yw (m)
L_inout*3/4 2*w_ch*1.25
L_inout*3/4 -2*w_ch*1.25
L_inout*3/4-5*w_ch*tan(ang_inout/2) -2*w_ch*1.25
Work Plane 1 (wp1)>Difference 1 (dif1)

1 In the Work Plane toolbar, click Booleans and Partitions and choose Difference.
2 Select the object arr1(1,3,1) only.
3 In the Settings window for Difference, locate the Difference section.
4 Click Clear Selection.
5 Select the objects arr1(1,1,1), arr1(1,2,1), and arr1(1,3,1) only.
6 Find the Objects to subtract subsection. Select the Activate Selection toggle button.
7 Select the objects c1 and pol1 only.
Work Plane 1 (wp1)>Polygon 2 (pol2)

1 In the Work Plane toolbar, click Polygon.
2 In the Settings window for Polygon, locate the Coordinates section.
xw (m) yw (m)
L_inout*3/4 2*w_ch*1.25
L_inout*3/4 -2*w_ch*1.25
L_inout*3/4+5*w_ch*tan(ang_inout/2) -2*w_ch*1.25

Work Plane 1 (wp1)>Difference 2 (dif2)
1 In the Work Plane toolbar, click Booleans and Partitions and choose Difference.
2 Select the objects arr1(1,1,2), arr1(1,2,2), and arr1(1,3,2) only.
3 In the Settings window for Difference, locate the Difference section.
4 Find the Objects to subtract subsection. Select the Activate Selection toggle button.
5 Select the object pol2 only.
Work Plane 1 (wp1)>Rotate 1 (rot1)

1 In the Work Plane toolbar, click Transforms and choose Rotate.
2 Select the object dif1 only.
3 In the Settings window for Rotate, locate the Rotation section.
4 In the Angle text field, type ang_inout.
5 Locate the Center of Rotation section. In the xw text field, type L_inout*3/4.
6 In the yw text field, type 2*w_ch*1.25.
Work Plane 1 (wp1)>Union 1 (uni1)

1 In the Work Plane toolbar, click Booleans and Partitions and choose Union.
2 Click in the Graphics window and then press Ctrl+A to select both objects.
3 In the Settings window for Union, click Build Selected.
Work Plane 1 (wp1)>Move 1 (mov1)

1 In the Work Plane toolbar, click Transforms and choose Move.
2 Select the object uni1 only.
3 In the Settings window for Move, locate the Displacement section.
4 In the xw text field, type -L_inout-w_ch*0.5.
5 In the yw text field, type -w_ch*2.5.

3 In the Width text field, type w_ch.

4 In the Height text field, type L_ch+10*w_ch.
5 Locate the Position section. In the xw text field, type -w_ch/2.
6 In the yw text field, type -5*w_ch.

3 Select the object r3 only.
5 In the xw size text field, type N_ch.
6 Locate the Displacement section. In the xw text field, type w_ch*2.

3 In the Width text field, type (N_ch-0.5)*(2*w_ch).
5 Locate the Position section. In the xw text field, type -w_ch/2.
6 In the yw text field, type -w_ch*5.

2 Select the object r4 only.
4 In the yw size text field, type 3.
5 Locate the Displacement section. In the yw text field, type w_ch*2.
Work Plane 1 (wp1)>Copy 1 (copy1)

1 In the Work Plane toolbar, click Transforms and choose Copy.
3 Select the objects arr3(1,1), arr3(1,2), and arr3(1,3) only.
4 In the Settings window for Copy, locate the Displacement section.

5 In the yw text field, type L_ch+5*w_ch.
Work Plane 1 (wp1)>Rotate 2 (rot2)

1 In the Work Plane toolbar, click Transforms and choose Rotate.
2 Select the object mov1 only.
4 In the Angle text field, type 180.
5 Locate the Center of Rotation section. In the xw text field, type (N_ch-1)*w_ch.
6 In the yw text field, type L_ch/2.
7 Locate the Input section. Select the Keep input objects check box.
Work Plane 1 (wp1)>Union 2 (uni2)

1 In the Work Plane toolbar, click Booleans and Partitions and choose Union.
2 Click the Select Box button in the Graphics toolbar.
3 Click in the Graphics window and then press Ctrl+A to select all objects.
Fillet Selection 1
1 In the Work Plane toolbar, click Selections and choose Box Selection.
2 In the Settings window for Box Selection, type Fillet Selection 1 in the Label text
field.
3 Locate the Geometric Entity Level section. From the Level list, choose Point.
4 Locate the Box Limits section. In the xw maximum text field, type 0.
5 In the yw minimum text field, type L_ch+w_ch*4.5.
Fillet Selection 2
1 In the Work Plane toolbar, click Selections and choose Box Selection.
2 In the Settings window for Box Selection, type Fillet Selection 2 in the Label text
field.
3 Locate the Geometric Entity Level section. From the Level list, choose Point.
4 Locate the Box Limits section. In the xw minimum text field, type 2*w_ch*(N_ch-1).
5 In the yw maximum text field, type -w_ch*4.5.

Work Plane 1 (wp1)>Fillet 1 (fil1)
1 In the Work Plane toolbar, click Fillet.
2 In the Settings window for Fillet, locate the Points section.
3 From the Vertices to fillet list, choose Fillet Selection 1.
4 Locate the Radius section. In the Radius text field, type w_ch/2.
Work Plane 1 (wp1)>Fillet 2 (fil2)

1 Right-click Component 1 (comp1)>Geometry 1>Work Plane 1 (wp1)>Plane Geometry>
Fillet 1 (fil1) and choose Duplicate.
2 In the Settings window for Fillet, locate the Points section.
3 Find the Vertices to fillet subsection. Select the Activate Selection toggle button.
4 From the Vertices to fillet list, choose Fillet Selection 2.
Extrude 1 (ext1)
1 In the Model Builder window, right-click Geometry 1 and choose Extrude.
2 In the Settings window for Extrude, locate the Distances section.
Distances (m)
h_a
Cylinder 1 (cyl1)
1 In the Geometry toolbar, click Cylinder.
2 In the Settings window for Cylinder, locate the Size and Shape section.
3 In the Radius text field, type R_in.
4 In the Height text field, type 3*h_a.
Rotate 1 (rot1)
1 In the Geometry toolbar, click Transforms and choose Rotate.
2 Select the object cyl1 only.
4 In the Angle text field, type ang_inout.

5 Locate the Point on Axis of Rotation section. In the x text field, type L_inout*3/4.
6 In the y text field, type 2*w_ch*1.25.
Move 1 (mov1)
1 In the Geometry toolbar, click Transforms and choose Move.
2 Select the object rot1 only.
3 In the Settings window for Move, locate the Displacement section.
4 In the x text field, type -L_inout-w_ch*0.5.
5 In the y text field, type -w_ch*2.5.
Rotate 2 (rot2)
1 In the Geometry toolbar, click Transforms and choose Rotate.
2 Select the object mov1 only.
3 In the Settings window for Rotate, locate the Input section.
4 Select the Keep input objects check box.
5 Locate the Rotation section. In the Angle text field, type 180.
6 Locate the Point on Axis of Rotation section. In the x text field, type (N_ch-1)*w_ch.
7 In the y text field, type L_ch/2.
Form Union (fin)

1 In the Model Builder window, click Form Union (fin).
2 In the Settings window for Form Union/Assembly, click Build Selected.
Inlet Manifold
1 In the Geometry toolbar, click Selections and choose Box Selection.
2 In the Settings window for Box Selection, type Inlet Manifold in the Label text field.
3 Locate the Box Limits section. In the x minimum text field, type -L_inout.
4 In the y maximum text field, type w_ch/2.
5 Locate the Output Entities section. From the Include entity if list, choose
Entity inside box.

Outlet Manifold
2 In the Settings window for Box Selection, type Outlet Manifold in the Label text field.
3 Locate the Box Limits section. In the x maximum text field, type N_ch*w_ch*2+
L_inout-w_ch.
4 In the y minimum text field, type L_ch-w_ch/2.
Entity inside box.
Channels Above Electrode Surface

2 In the Settings window for Box Selection, type Channels Above Electrode Surface
in the Label text field.
3 Locate the Box Limits section. In the y minimum text field, type -w_ch/2.
4 In the y maximum text field, type L_ch+w_ch/2.
Entity inside box.
Electrode Surface
2 In the Settings window for Box Selection, type Electrode Surface in the Label text
field.
3 Locate the Geometric Entity Level section. From the Level list, choose Boundary.
4 Locate the Box Limits section. In the y minimum text field, type -w_ch/2.
5 In the y maximum text field, type L_ch+w_ch/2.
6 In the z maximum text field, type h_a/2.
Entity inside box.
Inlet
2 In the Settings window for Box Selection, type Inlet in the Label text field.
4 Locate the Box Limits section. In the x maximum text field, type 0.
5 In the z minimum text field, type h_a*2.

Entity inside box.
Outlet
2 In the Settings window for Box Selection, locate the Geometric Entity Level section.
3 From the Level list, choose Boundary.
4 In the Label text field, type Outlet.
5 Locate the Box Limits section. In the x minimum text field, type N_ch*w_ch*2.
6 In the z minimum text field, type h_a*2.
Entity inside box.
Exterior Boundaries to Electrode Channels

1 In the Geometry toolbar, click Selections and choose Adjacent Selection.
2 In the Settings window for Adjacent Selection, type Exterior Boundaries to
Electrode Channels in the Label text field.
3 Locate the Input Entities section. Click Add.

4 In the Add dialog box, select Channels Above Electrode Surface in the Input selections list.
5 Click OK.
Exteror Boundaries to Inlet Manifold

2 In the Settings window for Adjacent Selection, type Exteror Boundaries to Inlet
Manifold in the Label text field.

4 In the Add dialog box, select Inlet Manifold in the Input selections list.
5 Click OK.
Inlets to Electrode Channels

1 In the Geometry toolbar, click Selections and choose Intersection Selection.
2 In the Settings window for Intersection Selection, type Inlets to Electrode
Channels in the Label text field.

5 In the Add dialog box, in the Selections to intersect list, choose
Exterior Boundaries to Electrode Channels and Exteror Boundaries to Inlet Manifold.
6 Click OK.
Manifolds
1 In the Geometry toolbar, click Selections and choose Union Selection.
2 In the Settings window for Union Selection, type Manifolds in the Label text field.
4 In the Add dialog box, in the Selections to add list, choose Inlet Manifold and
Outlet Manifold.
5 Click OK.
Exterior Boundaries to Manifolds

2 In the Settings window for Adjacent Selection, type Exterior Boundaries to
Manifolds in the Label text field.

4 In the Add dialog box, select Manifolds in the Input selections list.
5 Click OK.
Inlets and Outlets to Electrode Channels

2 In the Settings window for Intersection Selection, type Inlets and Outlets to
Electrode Channels in the Label text field.
Exterior Boundaries to Electrode Channels and Exterior Boundaries to Manifolds.
6 Click OK.
Electrode Channels and Manifolds

1 In the Geometry toolbar, click Selections and choose Union Selection.
2 In the Settings window for Union Selection, type Electrode Channels and Manifolds
in the Label text field.

4 In the Add dialog box, in the Selections to add list, choose Inlet Manifold, Outlet Manifold,
and Channels Above Electrode Surface.
5 Click OK.
Exteror Boundaries to Outlet Manifold

2 In the Settings window for Adjacent Selection, type Exteror Boundaries to Outlet
Manifold in the Label text field.
4 In the Add dialog box, select Outlet Manifold in the Input selections list.
5 Click OK.
Outlets from Electrode Channels

2 In the Settings window for Intersection Selection, type Outlets from Electrode
Channels in the Label text field.
Exterior Boundaries to Electrode Channels and Exteror Boundaries to Outlet Manifold.
6 Click OK.


Models - Cfd.pem Electrolyzer

Uploaded by

Copyright:

Available Formats

Models - Cfd.pem Electrolyzer

Uploaded by

Document Information

Original Description:

Original Title

Copyright

Available Formats

Share this document

Share or Embed Document

Sharing Options

Did you find this document useful?

Is this content inappropriate?

Copyright:

Available Formats

Models - Cfd.pem Electrolyzer

Uploaded by

Copyright:

Available Formats

Created in COMSOL Multiphysics 5.

Polymer Electrolyte Membrane Electrolyzer

2 | POLYMER ELECTROLYTE MEMBRANE ELECTROLYZER

Outlet (liquid water + oxygen gas)

Inlet (liquid water)

Figure 1: Geometry of the anode flow field.

3 | POLYMER ELECTROLYTE MEMBRANE ELECTROLYZER

Results and Discussion

Figure 2: Slice plot of the velocity magnitude at t = 10 s.

The Reynolds number Re is defined as

4 | POLYMER ELECTROLYTE MEMBRANE ELECTROLYZER

Figure 3: Gas volume fraction in the cell at t = 10 s.

5 | POLYMER ELECTROLYTE MEMBRANE ELECTROLYZER

Figure 4: Pressure drop in the cell at t = 0 s.

6 | POLYMER ELECTROLYTE MEMBRANE ELECTROLYZER

7 | POLYMER ELECTROLYTE MEMBRANE ELECTROLYZER

Figure 6: Individual channel velocities at various times.

Notes About the COMSOL Implementation

2. J. Nie, Y. Chen, S. Cohen, B. Carter, and R. Boehm, “Numerical and experimental

Application Library path: CFD_Module/Multiphase_Flow/pem_electrolyzer

8 | POLYMER ELECTROLYTE MEMBRANE ELECTROLYZER

1 In the Geometry toolbar, click Insert Sequence.

9 | POLYMER ELECTROLYTE MEMBRANE ELECTROLYZER

Name Expression Value Description

10 | POLYMER ELECTROLYTE MEMBRANE ELECTROLYZER

11 | POLYMER ELECTROLYTE MEMBRANE ELECTROLYZER

12 | POLYMER ELECTROLYTE MEMBRANE ELECTROLYZER

13 | POLYMER ELECTROLYTE MEMBRANE ELECTROLYZER

MIXTURE MODEL, LAMINAR FLOW (MM)

Inlet - Liquid Water

Inlet - Electrode Surface Oxygen Evolution

14 | POLYMER ELECTROLYTE MEMBRANE ELECTROLYZER

1 In the Model Builder window, under Component 1 (comp1) click Mesh 1.

15 | POLYMER ELECTROLYTE MEMBRANE ELECTROLYZER

16 | POLYMER ELECTROLYTE MEMBRANE ELECTROLYZER

17 | POLYMER ELECTROLYTE MEMBRANE ELECTROLYZER

18 | POLYMER ELECTROLYTE MEMBRANE ELECTROLYZER

Boundary Layer Properties

19 | POLYMER ELECTROLYTE MEMBRANE ELECTROLYZER

Boundary Layer Properties

20 | POLYMER ELECTROLYTE MEMBRANE ELECTROLYZER

Step 2: Time Dependent

21 | POLYMER ELECTROLYTE MEMBRANE ELECTROLYZER

22 | POLYMER ELECTROLYTE MEMBRANE ELECTROLYZER

23 | POLYMER ELECTROLYTE MEMBRANE ELECTROLYZER

Gas Volume Fraction

24 | POLYMER ELECTROLYTE MEMBRANE ELECTROLYZER

Velocity in Electrode Channels

25 | POLYMER ELECTROLYTE MEMBRANE ELECTROLYZER

Velocity in Electrode Channels

26 | POLYMER ELECTROLYTE MEMBRANE ELECTROLYZER

Name Expression Value Description

27 | POLYMER ELECTROLYTE MEMBRANE ELECTROLYZER

1 In the Model Builder window, under Component 1 (comp1)>Definitions click

28 | POLYMER ELECTROLYTE MEMBRANE ELECTROLYZER

1 In the Model Builder window, under Study 1>Solver Configurations right-click

Solution - Small Geometry

29 | POLYMER ELECTROLYTE MEMBRANE ELECTROLYZER