Assuming you have OpenFOAM version 6 with environment variables set, here is all you need to download, compile, and run sofcFoam.

obtain the code

cd <myChosenParentDirectory>
git clone git://git.code.sf.net/p/openfuelcell/git openfuelcell
cd openfuelcell

compile the model

cd src
./Allwmake
cd .. 

generate meshes

cd run/<caseDirectory> #coFlow, counterFlow, crossFlow, ...
make mesh

run model from the command line with delivered settings

make srun

generate VTK files for final output time

make view

The sofcFoam model is an OpenFOAM [1]application for the simulation of solid oxide fuel cells. It is a single-cell and multi-cell stack model being developed within the open source Multi-Scale Integrated Fuel Cell (MuSIC) project. This document describes how to obtain and use sofcFoam.

To make more sense of what follows, it may be useful to begin with a brief overview of the model. (A more detailed view will be provided later). The model simulates a four region fuel cell, as depicted schematically in Figure 1(a). Between two interconnects we find an air region, an electrolyte, and a fuel region. The model uses a computational domain for each of these regions, and one more global domain for the entire cell. Each domain supports its own fields. Pressure, momentum and species mass fractions, for example, are solved on the air and fuel domains, and temperature is solved on the global domain. Global and regional information is transferred back and forth via grid cell mappings that are established during mesh generation/splitting. Porous cathode and anode zones, shown in Figure 1(b), are incorporated into the air and fuel regions, respectively, using Darcy’s Law. Electrochemistry is assumed to occur on the electrode-electrolyte interfaces. The resulting electrochemical mass fluxes give rise to Dirichlet velocity conditions and Neumann mass-fraction conditions on the air and fuel boundaries interfacing the electrolyte.
After an initialization phase, the model enters an iteration loop as follows:

1. global temperature is mapped to the fluid regions, air and fuel
2. fluid densities are calculated
3. pressure and momentum are solved in fluid domains
4. mass diffusivity in fluids is calculated
5. mass fractions of fluid species are solved
6. electrochemistry is calculated
7. regional velocity and thermophysical data fields are mapped to global mesh
8. energy equation is solved for global temperature
9. Steps (1) to (8) are repeated until convergence

OpenFoam
The openFuelCell program is being developed with the open source C++ objected oriented software OpenFOAM, therefore it is necessary to have a working installation of OpenFOAM.  At time of writing, OpenFuellCell has been tested and is known to work with the repository release of 2.1.x [2].

Some experience running OpenFOAM applications, as might be obtained from OpenFoam tutorials, will be helpful. Instructions for building OpenFOAM come with the download, and further assistance can be obtained from the discussion forum [3]. The WIKKI pages [4]are also useful for supplementary information.

Git
The sofcFoam code is maintained in a Git  [5] version control system repository on SourceForge [6].  The Git command line client tool is git.  Graphical tools are also available, e.g. , SmartGit [7].

A tar file containing the latest version of the OpenFuelCell code can be downloaded here from SourceForge [8]. 

Alternatively, users can check out the latest version of the code from the SourceForge git repository using the command...

git clone git://git.code.sf.net/p/openfuelcell/git openfuelcell

In your chosen parent directory for the sofcFoam model, (e.g. your OpenFoam work space $WM_PROJECT_USER_DIR/applications/), check out the master branch, from the git repository [9] using your favourite graphical git tool, or the command

git clone git://git.code.sf.net/p/openfuelcell/git openfuelcell

This creates directory openfuelcell/ in the current working directory.

src
To compile the library and application source code, go to openfuelcell/src/ directory and run the Allwmake script. Type ./Allwmake at the prompt. This should generate shared object library libsofcFoam.so in the $FOAM_USER_LIBBIN directory and application executable sofcFoam in the $FOAM_USER_APPBIN directory. A lnInclude/ directory, containing links to all of the libSrc class files, will appear in the libSrc/ directory.

cases
As can be seen in Figure 2(b), a case directory contains only one polyMesh/ directory immediately after checkout, and it contains only the dictionary file blockMeshDict. This dictionary, together with the setSet batch command files in the <case>/config/ directory, describes the global and regional meshes. After the global mesh is made by the OpenFOAM utility blockMesh, the utility splitMeshRegions generates the required regional meshes and map files. For more information on the blockMesh, setSet, and setsToZones utilities, see Chapter 5 “Mesh generation and conversion [10]” and Section 3.6 “Standard utilities [11]” in the OpenFOAM User Guide [12].

Making the global and regional meshes is handled in sofcFoam by the Makefile in the case directory. See, for example, run/coFlow/Makefile. The command

make mesh

issued from the case directory, will generate the global mesh and the region meshes. During model execution, various material property and other field values will be mapped from the region meshes to the global mesh. Cells that began life labeled as a fluid in the global mesh may have become a solid, and some of these may have boundary faces on the original fluid inlet or outlet patches. Accordingly, the fluid inlet and outlet patches may need to be redefined for the new reality. The redefinitions are specified by the make.face[Air|Fuel|Set] files in the config directory. See Appendix A for a description of the steps required to specify a new geometry.

Running on a sequential machine:

With the application compiled, the command:

make srun

will run the executable from the command line, using the available case data. The model can also be run by typing the executable name, and the output directed to Standard Out can be redirected to a file:

sofcFoam | tee log.run

Instead of running the model from the command line, a runscript is available to submit a job to a queue. The script usage line may need editing for your queuing system.

After the model has run to completion, VTK files for visualization, e.g. with paraview, can be prepared easily using the Makefile. Typing:

make view

will generate VTK files for the last output step, whereas:

make viewAll

will generate VTK files for all output directories.

Running on a parallel machine:

Case quickTestPar is essentially case quickTest, modified for parallel decomposition. The makefile, however, is different from the serial file. 

Edit system/decomposeParDict

Export environment variable NPROCS to desired number of processors (e.g. in csh: setenv NPROCS n)

make mesh

Builds the mesh

make parprep

Decompose the mesh and initial field data

make run

Run the decomposed case in parallel

make reconstruct

Reconstruct the mesh and data

make view

Generate the VTK files for viewing.

To start over type:

./Allclean

Here, the field files in the region subdirectories 0/air and 0/fuel have counterparts on the global mesh, i.e., in the top level of the 0/ directory. These files in the top level 0/ directory will be subset during parallel mesh splitting, and the subset fields will be moved to the processor/0/air and processor/0/fuel directories.

In the out−of−box case configuration, air species are O2 and N2, while fuel species are H2 and H2O. If species Sp are added or changed, it is necessary to have diffSp and YSp files for the new species, in accordance with sections 11.8.4 and 11.8.5 of the "gettingStarted_Cell268_OF21x.pdf" document in the docs/ directory. In order to reconstruct and visualize the fluid region fields, the fluid region polyMesh subdirectories must be present in the constant directory. Furthermore, the *ProcAddressing files usually generated by decomposePar −region <name> are required in the processor polyMesh directories.

We also need the *RegionAddressing files in our processor polyMesh directories to enable mapping between global and region meshes.

Accomplishing all of the above requires two decompositions for parallelism, one for the region meshes only and one for the entire mesh. The second decomposition can’t proceed in the presence of the processor* directories created from the first decomposition, so we rename them to keep their *ProcAddressing files for later use. There are also two mesh splittings, the second of them in parallel.

We proceed with steps 1 to 11 below. Note that Step 2 is controlled by the makefile:

$ make mesh | tee log.mesh

Steps 3 through 11 are controlled by parprep.csh, either directly or through the Makefile:

$ ./parprep.csh or $ make parpre
Step 10 is controlled by cleanup.csh, which contains lists noGlobal, noAir and noFuel, of files that are to be removed from, respectively, the processor*/0, processor*/0/air and processor*/0/fuel directories. These lists may need to be edited for additions or changes, particularly with respect to species. cleanup.csh is executed within parprep.csh.

1. Edit system/decomposeParDict for the desired mesh decomposition and copy to system/air and system/fuel.
   Set environment variable NPROCS to number of processors in decomposeParDict
2. Build the mesh and split off the regions
3. Decompose the fluid regions for parallel run
4. Rename the processor* directories as proc_*
5. Decompose the global mesh in parallel
6. Split the global mesh in parallel
7. Install newly created region/polyMesh dirs and field files to processors
8. Mesh splitting introduces interface patches named <regionA>_to_<regionB>.
   Their field subdirectories need to be edited. A sed script exists for this.
9. Sets and zones need to be recreated for the processor* meshes
10. Directories are cleaned up by removing extraneous files
11. The *ProcAddressing files are copied from proc_* to corresponding dirctories in processor*

To run the code:

mpirun −mca mpi_warn_on_fork 0 −np $NPROCS sofcFoam −parallel

The script

prun

will issue the above command, and

make runwill launch the prun script.

To reconstruct the parallel run:

reconstructPar
reconstructPar −region air
reconstructPar −region fuel

or simply

make reconstruct

To make the VTK files

foamToVTK −latestTime
foamToVTK −latestTime −region air
foamToVTK −latestTime −region fuel
foamToVTK −latestTime −region electrolyte
foamToVTK −latestTime −region interconnect0
foamToVTK −latestTime −region interconnect1

or simply

make view

In this section, we will follow the main function, in the file sofcFoam.C, through the model execution. Like all OpenFOAM applications, the model begins by including the OpenFOAM src files setRootCase.H to check the case path, and createTime.H to read the system/controlDict file and instantiate the Time object runTime. This is followed by the creation of meshes, reading of properties, and creation of fields for the global cell mesh and the region meshes.

Figure A1 shows the proposed geometry we intend to model. The associated dimensions of the components are given in Table A1. The vertical structure can be captured by seven blocks, as shown in Figure A2, (the block containing the electrolyte is too thin to be discernible). The blocks containing the air and fuel channels can then be split horizontally to separate the channels from the ribs, the latter being part of the interconnects.

Figure A1

Figure A1. A fuel cell with one air channel and one fuel channel. Left panel shows air (blue) and fuel (purple) inlets, interconnects (grey) and electrode sides. Centre panel shows air (blue) and fuel (purple) volume regions, each comprised of both a channel and a porous electrode zone. Right panel shows lower (blue) and upper (red) interconnect regions.

We begin with a blockMeshDict dictionary that will create a parent mesh consisting of the seven vertical blocks (Figure A2), which for convenience, going from bottom to top, we refer to as interconnect0, air, cathode, electrolyte, anode, fuel, and interconnect1. Although the geometry shows symmetry about the y = 2 plane, we construct the entire domain for illustrative purposes. Here is the list of points for the blockMeshDict file:

blockMeshDict
convertToMeters 0.001;
vertices
(
// ... From Bottom To Top
// Interconnect0
( 0 0 0) // 0
(50 0 0) // 1
(50 4 0) // 2
( 0 4 0) // 3
// Interconnect0_to_Air
( 0 0 3.5) // 4
(50 0 3.5) // 5
(50 4 3.5) // 6
( 0 4 3.5) // 7
// Air_to_cathode
( 0 0 5.0) // 8
(50 0 5.0) // 9
(50 4 5.0) //10
( 0 4 5.0) //11
// cathode_to_Electrolyte
( 0 0 5.29) //12
(50 0 5.29) //13
(50 4 5.29) //14
( 0 4 5.29) //15
// Electrolyte_to_anode
( 0 0 5.3) //16
(50 0 5.3) //17
(50 4 5.3) //18
( 0 4 5.3) //19
// anode_to_Fuel
( 0 0 6.3) //20
(50 0 6.3) //21
(50 4 6.3) //22
( 0 4 6.3) //23
// fuel_to_Interconnect1
( 0 0 7.8) //24
(50 0 7.8) //25
(50 4 7.8) //26
( 0 4 7.8) //27
// Interconnect1
( 0 0 11.3) //28
(50 0 11.3) //29
(50 4 11.3) //30
( 0 4 11.3) //31
);In the vertices section above, each set of four vertices defines a horizontal rectangle representing an interface between the above mentioned blocks (and including the top and bottom surfaces). As can be readily seen, the vertices of a rectangle are arranged so that a traversal from one to the next takes one anticlockwise around the rectangle, starting from x=0. Note that the coordinates are scaled by 0.001 metres, so the maximum x-coordinate, for example, is 50 mm. The vertices are numbered by their index in the list, beginning at index 0. Figure A3 shows the location of some of these points on the geometry.

Figure A3. Location on geometry of selected vertices, as numbered by the blockMeshDict file.

In the blocks section below, each hexahedral block is defined by two successive sets of four vertices, i.e. the corner vertices of the block. The air block, eg, is defined by: hex (4 5 6 7 8 9 10 11). The number of cells in each coordinate direction and the grading of the mesh are also prescribed here.

blocks
(
// Interconnect0
hex (0 1 2 3 4 5 6 7) (25 8 7) simpleGrading (1 1 1)
// air
hex (4 5 6 7 8 9 10 11) (25 8 3) simpleGrading (1 1 1)
// cathode
hex (8 9 10 11 12 13 14 15) (25 8 1) simpleGrading (1 1 1)
// electrolyte
hex (12 13 14 15 16 17 18 19) (25 8 1) simpleGrading (1 1 1)
// anode
hex (16 17 18 19 20 21 22 23) (25 8 2) simpleGrading (1 1 1)
// fuel
hex (20 21 22 23 24 25 26 27) (25 8 3) simpleGrading (1 1 1)
// Interconnect1
hex (24 25 26 27 28 29 30 31) (25 8 7) simpleGrading (1 1 1)
);We have no need to define any edges.
edges
(
);
A patch consists of one or more outer boundaries of the blocks. These boundaries (rectangles in our case) are described by their corner vertices, arranged so that a traversal from one to the next takes one round the rectangle anticlockwise about the outward normal.

patches
(
// ... From Bottom to Top
// Interconnect0
patch interconnect0Bottom
(
(0 3 2 1)
)
patch interconnect0Sides
(
(0 1 5 4)
(3 7 6 2)
(0 4 7 3)
(1 2 6 5)
)
// Air
patch airInlet
(
(4 8 11 7)
)
patch airOutlet
(
(5 6 10 9)
)
patch airSides
(
(4 5 9 8)
(7 11 10 6)
)
// Cathode
patch cathodeSides
(
(8 9 13 12)
(11 15 14 10)
(8 12 15 11)
(9 10 14 13)
)
// Electrolyte
patch electrolyteSides
(
(12 13 17 16)
(15 19 18 14)
(12 16 19 15)
(13 14 18 17)
)
// Anode
patch anodeSides
(
(16 17 21 20)
(19 23 22 18)
(16 20 23 19)
(17 18 22 21)
)
// Fuel
patch fuelInlet
(
(20 24 27 23)
)
patch fuelOutlet
(
(21 22 26 25)
)
patch fuelSides
(
(20 21 25 24)
(23 27 26 22)
)
// interconnect1
patch interconnect1Sides
(
(24 28 31 27)
(25 26 30 29)
(24 25 29 28)
(27 31 30 26)
)
patch interconnect1Top
(
(28 29 30 31)
)
);
mergePatchPairs
(
);
// **********************************

For more description of the blockMeshDict dictionary and the blockMesh utility, see section 5.3, Mesh generation with the blockMesh utility, in the OpenFoam User Guide [13].

We must now define the cellSets that will make up the cells of our five regions: interconnect0, air, electrolyte, fuel and interconnect1. Using the cellSets, a mesh will be generated for each region. Note that the cathode and anode blocks will become porousZones within the air and fuel regions, respectively. A portion of the air block contains two ribs that must become part of the interconnect0 region, and similarly two ribs contained in the fuel block must become part of the interconnect1 region. Cells in the electrolyte block will form the electrolyte region, and cells in the interconnect blocks will become part of the interconnect regions.

The cellSets for the regions are specified in config/make.setSet. Here each cellSet is defined by the diagonally opposite corners of a box bounded by coordinate planes.

The first set specified is the cellSet interconnect0. The specification begins with the cells in the interconnect0 block, which consists of all the cells below z=3.5mm (note that the coordinates are given in metres). Then the cells of the ribs are added. One of these extends from y=0 mm to y=1 mm, and the other from y=3 mm to y=4 mm. Both extend the full length of 50 mm in x, and in height from z=3.5 mm to z=5 mm.

The specification for the air cellSet begins with the cathode block and adds the channel, which extends the full length of 50 mm in x, from y=1 mm to y=3 mm in width, and from y=3.5 mm to y=5 mm in height. The remaining sets are similarly specified.

make.setSet

cellSet interconnect0 new boxToCell (0 0.0e-3 0.0e-3) (50.0e-3 4.0e-3 3.5e-3)
cellSet interconnect0 add boxToCell (0 0.0e-3 3.5e-3) (50.0e-3 1.0e-3 5.0e-3)
cellSet interconnect0 add boxToCell (0 3.0e-3 3.5e-3) (50.0e-3 4.0e-3 5.0e-3)
cellSet air new boxToCell (0 0.0e-3 5.0e-3) (50.0e-3 4.0e-3 5.29e-3)
cellSet air add boxToCell (0 1.0e-3 3.5e-3) (50.0e-3 3.0e-3 5.0e-3)
cellSet electrolyte new boxToCell (0 0 5.29e-3) (50.0e-3 4.0e-3 5.3e-3)
cellSet fuel new boxToCell (0 0.0e-3 5.3e-3) (50.0e-3 4.0e-3 6.3e-3)
cellSet fuel add boxToCell (0 1.0e-3 6.3e-3) (50.0e-3 3.0e-3 7.8e-3)
cellSet interconnect1 new boxToCell (0 0.0e-3 7.8e-3) (50.0e-3 4.0e-3 11.3e-3)
cellSet interconnect1 add boxToCell (0 0.0e-3 6.3e-3) (50.0e-3 1.0e-3 7.8e-3)
cellSet interconnect1 add boxToCell (0 3.0e-3 6.3e-3) (50.0e-3 4.0e-3 7.8e-3)

The air and fuel regions are each given a porous zone within the fluid zone, as specified in config/make.setAir and config/make.set fuel:

make.setAir

cellSet air new boxToCell (0 0.0e-3 5.0e-3) (50.0e-3 4.0e-3 5.29e-3)
cellSet air add boxToCell (0 1.0e-3 3.5e-3) (50.0e-3 3.0e-3 5.0e-3)
cellSet cathode new boxToCell (0 0 5.0e-3) (40.0e-3 4.0e-3 5.29e-3)

make.setFuel

cellSet fuel new boxToCell (0 0.0e-3 5.3e-3) (50.0e-3 4.0e-3 6.3e-3)
cellSet fuel add boxToCell (0 1.0e-3 6.3e-3) (50.0e-3 3.0e-3 7.8e-3)
cellSet anode new boxToCell (0 0 5.3e-3) (50.0e-3 4.0e-3 6.3e-3)

Clearly, the fluid inlet and outlet patches on the global mesh are incorrect, since their original definitions include faces that are really part of the interconnect ribs. The correction proceeds in three steps. First, faceSets for all of the existing patches of the blockMesh are created using the patchToFace action of the faceSet utility, as specified by the config/make.faceSet file:

faceSet interconnect0Sides new patchToFace interconnect0Sides all
faceSet interconnect0Bottom new patchToFace interconnect0Bottom all
faceSet interconnect1Sides new patchToFace interconnect1sides all
faceSet interconnect1Top new patchToFace interconnect1Top all
faceSet electrolyteSides new patchToFace electrolyteSides all
faceSet cathodeSides new patchToFace cathodeSides all
faceSet airSides new patchToFace airSides all
faceSet airInlet new patchToFace airInlet all
faceSet airOutlet new patchToFace airOutlet all
faceSet anodeSides new patchToFace anodeSides all
faceSet fuelSides new patchToFace fuelSides all
faceSet fuelInlet new patchToFace fuelInlet all
faceSet fuelOutlet new patchToFace fuelOutlet all
faceSet interconnect0Sides add patchToFace airInlet all
faceSet interconnect0Sides add patchToFace airOutlet all
faceSet interconnect0Sides add patchToFace airSides all
faceSet interconnect1Sides add patchToFace fuelInlet all
faceSet interconnect1Sides add patchToFace fuelOutlet all
faceSet interconnect1Sides add patchToFace fuelSides all
faceSet airSides clear
faceSet airInlet clear
faceSet airOutlet clear
faceSet fuelSides clear
faceSet fuelInlet clear
faceSet fuelOutlet clearNote that the make.faceset file also specifies some manipulations, adding faceSets airInlet, airOutlet, and airSides to the faceSet interconnect0, and similary on the fuel side. After being added, they are subsequently cleared. Next, the inlet and outlet faceSets are corrected using new specifications in config/make.faceAir and config/make.faceFuel:

make.faceAir

faceSet airInlet new boxToFace (-1e-6 1.0e-3 3.5e-3) (1e-6 3.0e-3 5.0e-3)
faceSet airOutlet new boxToFace (39.999e-3 1.0e-3 3.5e-3) (40.001e-3 3.0e-3 5.0e-3)
faceSet interconnect0Sides delete faceToFace airInlet all
faceSet interconnect0Sides delete faceToFace airOutlet all

make.faceFuel

faceSet fuelInlet new boxToFace (-1e-6 1e-3 6.3e-3) (1e-6 3.0e-3 7.8e-3)
faceSet fuelOutlet new boxToFace (39.999e-3 1e-3 6.3e-3) (40.001e-3 3.0e-3 7.8e-3)
faceSet interconnect1Sides delete faceToFace fuelInlet all
faceSet interconnect1Sides delete faceToFace fuelOutlet all

The new inlet and outlet patches are defined by a bounding box for the new patch. Here the new airInlet, eg, is normal to the x-direction and is bounded by a box which is shallow in x, extending 1e-6 m in front of and behind the prescribed x-coordinate location. The lateral extents of the box in the other two directions correspond to the lateral extent of the inlet in those directions. Faces with face centre within the box will be selected, so the box must not extend to the adjacent grid cell. The fuelInlet and the two outlets are similarly defined. The new inlets and outlets are then removed from the interconnect faceSets.

Finally, the faceSets are used to create new patches using the createPatch utility, which is controlled by the system/createPatchDict file. Here is an excerpt for the airInlet patch:

patchInfo
(
{
name airInlet;
// Type of new patch
dictionary
{
type patch;
}
constructFrom set;
patches ();
set airInlet;
}
. . .
);

We will find the following entry (with additional face numbering information) for the airInlet in the mesh boundary file.
airInlet
{
type patch;
}
The remaining patches are formed in the same way. The complete patch list is:
interconnect0Bottom
interconnect0Sides
airInlet
airOutlet
cathodeSides
electrolyteSides
anodeSides
fuelInlet
fuelOutlet
interconnect1Sides
interconnect1Top