7. Introduction to Microscale Simulations¶
Microstructure simulations have been introduced in OpenFCST release v0.3. OpenFCST provides tools to generate meshes from image stacks and then perform diffusion, diffusion with reaction boundary condition and electron transport simulations on the generated meshes. OpenFCST uses a VTK mesh generator using the TVTK library [1]. This example illustrates the different steps from generating a mesh from images to setting up and running simulations.
The steps covered in this tutorial are:
- Generating a VTK mesh using PythonFCST
- Diffusion based simulation in OpenFCST
- Electron transport simulation in OpenFCST
7.1. Generating a VTK mesh using PythonFCST¶
In this example, we will demonstrate how to generate an Unstructured Grid mesh in the Legacy VTK format using the PythonFCST package (which are Python libraries to extend functionality of OpenFCST). This utility to generate a VTK mesh is beneficial when a mesh needs to be created from a stack of tiff images (such as those obtained from FIB-SEM, TEM, etc.) or from a 3D numpy array in Python. There are two approaches that can be used which are described below
7.1.1. Using the Python executable writeVTK.py¶
PythonFCST comes equipped with a Python executable writeVTK.py which reads in a stack of tiff images and generates a Legacy Unstructured Grid file using the pixel values in the images as the material ids for the mesh. Before running writeVTK.py it is important that OpenFCST is installed on the system and the OpenFCST environment variables are set using the fcst_env.sh file in the FCST Install Directory (${FCST_DIR}). To set the environment variables type the following command in the Terminal (or Konsole) from ${FCST_DIR}:
$ . fcst_env.sh
Once the environment variables have been set the writeVTK.py should be available as a command line utility. To ensure the environment variable has been set, the following command can be used:
$ writeVTK.py -h
Usage: writeVTK.py [options]
Options:
-h, --help show this help message and exit
-f FILE, --file=FILE relative path of the images to be converted
-b BASENAME, --basename=BASENAME
basename of the image files. Default: 'Slice_' would
mean images with name Slice_001 and so on
-n NUM Number of images to be read. Default: 200
-o OUTPUT Output file name. Default: test.vtk
-m MATERIAL, --material=MATERIAL
material id of the phase to be meshed. Default: 0
-v VOXEL, --voxel=VOXEL
voxel size as x y z
-e EXTENSION, --extension=EXTENSION
extension for the input files. Default is .tiff
The above command generates the help documentation for the writeVTK.py. In order to generate the mesh, the relative folder location from the current working directory for the images needs to be specified with the -f identifier.
The basename of the images refers to the name of the image without the numerical part and extension. For example, for images named ‘Slice_001.tiff’ to ‘Slice_100.tiff’ the basename would be 'Slice_'. The number of images to be read in must be specified with the -n identifier.
The output file name can be specified with the -o identifier. The voxel size which is the dimension in cm per pixel in the X,Y and Z direction is specified with the -v identifier.
By default the extension for the images is ‘.tiff’. If the images have a different extension then they must be specified with the -e identifier.
Another important option is the material id which can be specified with the -m option. The material id represents the pixel value for which the mesh should be generated. The default material id is 0 which corresponds to the black color in a grey-scale image.
For our example, we will generate a mesh from 10 images in the folder ${FCST_DIR}/examples/microscale/3D_Diffusion/images while the current working directory is ${FCST_DIR}/examples/microscale/3D_Diffusion, the following command is used:
$ writeVTK.py -f images/ -b Slice_ -n 10 -o my_test.vtk -m 0 -v 1e-6 1e-6 1e-6 -e .tiff
==================================================
==================================================
= Points and cells Created
= Writing a VTK file
==================================================
Calculating Fitted Sphere Size Distribution
**************************************************
= Processing dataset
--------------------------------------------------
==================================================
= Calculate PSD
= - Pore radius 1e-06
= - Pore radius 1.1643739803e-06
= - Pore radius 1.32874796059e-06
= - Pore radius 1.49312194089e-06
= - Pore radius 1.65749592118e-06
= - Pore radius 1.82186990148e-06
= - Pore radius 1.98624388177e-06
= - Pore radius 2.15061786207e-06
= - Pore radius 2.31499184237e-06
= - Pore radius 2.47936582266e-06
= - Pore radius 2.64373980296e-06
= - Pore radius 2.80811378325e-06
= - Pore radius 2.97248776355e-06
= - Pore radius 3.13686174384e-06
= - Pore radius 3.30123572414e-06
= - Pore radius 3.46560970443e-06
= - Pore radius 3.62998368473e-06
= - Pore radius 3.79435766503e-06
= - Pore radius 3.95873164532e-06
= - Pore radius 4.12310562562e-06
= The file has been successfully written!
--------------------------------------------------
Phase with material id: 0 is represented by material id 0 in the mesh file: my_test.vtk
time elapsed is 17.7439901829 seconds
If the command runs successfully, the above output will be displayed and the mesh file named my_test.vtk will be written to the ${FCST_DIR}/examples/microscale/3D_Diffusion (which was the current working directory). In addition to generating the mesh, the code also calculates the local radius based on a sphere fitting algorithm [2] of the phase and stores it in the generated mesh which will be later used in OpenFCST to account for the local Knudsen effects in the diffusion simulations.
7.1.2. Using the MESH module in PythonFCST¶
For people who are familiar with Python, using the PythonFCST.mesh module might be another option to generate the mesh. The mesh module contains the Grid Generator class which generates the VTK mesh. The arguments to be supplied are a 3D numpy array (argument name: image) and the output filename (argument name: filename). This provides a bit more flexibility to use any 3D numpy array which might have been generated from any source rather than being limited to a stack of tiff images.
For the purpose of this example, we would deal with the same set of images as above. Let’s use the PythonFCST.mesh module to generate the mesh. The code snippet below describes how to use it for the example,
#Initial imports for libraries and modules
import PythonFCST.mesh.PhaseGenerator as grid
from PIL import Image
import os
import numpy as np
#Defining variables to read in images; Omit this section if you already have a 3D array
foldername = "images"
basename = "Slice_"
num_images = 10
output = "my_test.vtk"
name=os.getcwd()+'/'+foldername+'/'+basename+'%.3d'
extension='.tiff'
tmp = []
#Looping over images and reading them
for i in range(num_images):
filename=name%(i+1)+extension
tmp.append(np.asarray(Image.open(filename),dtype=np.float))
#Generating the 3D numpy array
image=np.dstack(tmp)
image[image>0]=20 #changing all non-zero material ids to 20; assuming that pixel value of 0 was the one of interest
#Generating the VTK mesh using the Mesh module in PythonFCST
o=grid(image,output,material=0)
o.write()
The generated mesh file: my_test.vtk has the material id 0 from the grey scale images extracted.
Mesh of the pores represented by the black color (pixel value 0) in images
Representation of the boundary ids [1-6] for the generated mesh
The mesh above shows the Knudsen radius of the pores computed and stored in the mesh. The boundary ids in the X, Y and Z direction are 1-2, 3-4 and 5-6 respectively. It is important to know these boundary ids because they will be used in the simulations.
Note
VTK files invert the X and Y axis from the images. Therefore, if the image had dimensions 20 x 30 in X and Y direction respectively, then the mesh will have dimensions 30 x 20 in the X and Y direction. By default, OpenFCST corrects for this so that the voxel sizes entered by user are in the true X and Y directions. Also, the boundary ids are corrected so that 1-2 are boundary ids in the true X direction and in the Y direction for the mesh as seen above.
The mesh can be further pre-processed using any of the filters in Paraview or Mayavi and could be written to a Legacy VTK file. This would still be compatible with the OpenFCST application.
7.2. Diffusion based simulation in OpenFCST¶
7.2.1. Introduction¶
Diffusion based simulations i.e, molecular diffusion, Knudsen diffusion and diffusion with reaction, are implemented in OpenFCST in the application AppDiffusion. The application is used to simulate gas diffusion in a defined phase, i.e., pores in the microstructures, and can be used to compute the effective diffusivity of the porous media. In addition, the application can also be used to simulate oxygen diffusion with a reaction boundary condition and predict the electrochemical performance in the microstructure assuming that the overpotential losses are negligible.
7.2.2. Governing equation¶
The governing equation for the application is based on Fick’s Law and given by,
where 
is the oxygen diffusion coefficient, 
is the total gas concentration (assumed constant), and 
is the molar fraction of oxygen which is also the solution variable for this application.
If Knudsen effects are not taken into account then is the binary diffusion coefficient of oxygen in air.
If Knudsen effects are to be taken into account then is modified to account for the Knudsen effects using the Bosanquet equation,
where 
is the molecular diffusion coefficient of the oxygen gas mixture (which can be supplied by user or calculated using the Chapman-Enskog theory) and 
is the Knudsen diffusion coefficient computed as follows,
where 
is the local pore radius obtained using the pore size distribution algorithm (shown above), 
is the universal gas constant, 
is the temperature of the gas and 
is the molecular weight of oxygen.
Note
is computed internally during execution based on the KnudsenRadius field supplied in the VTK mesh. Currently, only VTK meshes generated using the OpenFCST meshers are supported for accounting for local Knudsen effects in microstructures.
There are two sets of boundary conditions that can be used with AppDiffusion. The first set of boundary conditions are for the case of pure diffusion and given by,
where 
is the inlet plane and 
is the outlet plane opposite to the inlet plane. The governing equation with these BCs are linear and therefore, can be solved using direct solvers such as UMFPACK, MUMPS or iterative solvers such as CG, GMRES.
The second of boundary conditions are for the case when reaction source terms are added to the oxygen diffusion. For this case the boundary conditions are as follows,
where 
is the current density, 
is a ratio of the area of active platinum surface to the interfacial area between solid and pore and 
denotes the solid-pore interface.
For more details on the calculation of the source term please see [2]. The boundary conditions and the governing equations form a non-linear system of equations which are linearized using Picard iterations and inner loops are solved using a linear solver such as MUMPS, UMFPACK or CG.
7.2.3. Setting up a diffusion simulation (with & without Knudsen effects)¶
Once we have the mesh generated using either the above method or any of the other methods, it is time to setup the simulation in OpenFCST. OpenFCST has the Fick’s binary diffusion model implemented as the application diffusion. This can be used for simulating the flow through one of the phases in the 3D mesh that has been generated earlier. Before, running the simulation it is necessary that the environment variables are updated for the OpenFCST settings as described earlier. Once the environment is setup, we need to setup the parameter files for running the simulation.
As a first example, we setup a simple diffusion simulation without Knudsen effects. This example is contained in the folder ${FCST_DIR}/examples/microscale/3D_Diffusion For this the first file is the main.prm which is the argument file supplied to the OpenFCST executable. The main.prm file should look like this:
###################################################################### # $Id: main.prm 2011-02-09 secanell $ # # # # This file is used to simulate app_diffusion and obtain a # # a concentration profile # # # # # # Copyright (C) 2011 by Marc Secanell # # # ###################################################################### ###################################################################### subsection Simulator set simulator name = diffusion set simulator parameter file name = data.prm set nonlinear solver name = None set Analysis type = Analysis end ######################################################################
The main.prm sets up the overall simulation. The parameters that need to be specified for this tutorial are the simulator name, to diffusion, and the type of solver, to Linear. Since the Fick’s binary diffusion model results in a linear system of equations it does not require the use of non-linear solvers. In addition we specify the simulator parameter file to data.prm which specifies the necesssary parameters required to carry out the simulation. The data.prm looks like this:
######################################################################
# $Id: data.prm 2011-02-09 secanell $ #
# #
# This file is used to simulate app_diffusion and obtain a #
# a concentration profile. This file supplies the necessary #
# parameters for the simulation. #
# #
# Copyright (C) 2011 by Marc Secanell #
# #
######################################################################
###############
subsection Grid generation
set Type of mesh = GridExternal # Cathode | CathodeMPL | GridExternal
set File name = my_test.vtk
set File type = vtk
set Initial refinement = 0
set Sort Cuthill-McKee = false
end
###############
subsection Adaptive refinement
set Refinement = global #global | adaptive
set Number of Refinements = 2
set Output intermediate solutions = true
set Output final solution = true
set Output intermediate responses = true
set Use nonlinear solver for linear problem = false
end
###############
subsection System management
set Number of solution variables = 1
subsection Solution variables
set Solution variable 1 = oxygen_molar_fraction
end
subsection Equations
set Equation 1 = Ficks Transport Equation - oxygen
end
end
###############
subsection Equations
subsection Ficks Transport Equation - oxygen
subsection Initial data
set oxygen_molar_fraction = 0: 0.1 # where 0 indicates the material_id setup in the grid and 0.1 is the oxygen molar fraction
end
subsection Boundary data
set oxygen_molar_fraction = 3: 0.4, 4:0.01 # where 3 & 4 denote the boundary ids and 0.4 and 0.01 are the oxygen molar fraction
end
end
end
###############
subsection Discretization
set Element = FESystem[FE_Q(1)] #FESystem[FE_Q(3)-FE_Q(1)^2] #FESystem[FE_Q(1)^3] #System of three fe
subsection Matrix
set Quadrature cell = -1
set Quadrature face = -1
end
subsection Residual
set Quadrature cell = -1
set Quadrature face = -1
end
end
###############
subsection Fuel cell data
####
subsection Operating conditions
set Temperature cell [K] = 353
set Cathode pressure [Pa] = 101265 # (1 atm)
set Cathode relative humidity = 0.6
end
####
####
subsection Cathode gas diffusion layer
set Gas diffusion layer type = DummyGDL #[ DesignFibrousGDL | DummyGDL | SGL24BA ]
set Material id = 0
####
subsection DummyGDL
set Oxygen diffusion coefficient, [cm^2/s] = 0.22
end
####
end
####
end
###############
subsection Output Variables
set Compute boundary responses = true
set num_output_vars = 1
set Output_var_0 = oxygen_molar_fraction
end
###############
subsection Output
subsection Data
set Output format = vtu
set Print solution = true
end
subsection Grid
set Format = eps
end
end
################################
The data.prm comprises of several important parameters which can be modified to affect the simulation. The first section is the Grid generation where the Type of mesh needs to be specified as GridExternal because we are supplying a user-defined mesh and also, the File name & File type are specified based on the mesh generated earlier.
The next section is the Adaptive refinement where the option can be set for number of refinements in the solution as well as the various file output options can be changed.
The next section is System management where information with respect to the equations and solution variables is provided. For the Fick’s diffusion model, we will just have one solution variable. However, the gas to be simulated can be modified by changing the name of the solution variable and the equation. For list of gases that can be simulated refer to PureGas.cc. For our example, we use a oxygen diffusion in nitrogen (default) which is specified as the solution variable and equation.
Equations section allows to input the initial and boundary conditions for the model. Although for the diffusion model the initial data is not necessary but it could be supplied just to ensure that the mesh is being read in correctly and appropriate phase is being solved for. For this example, we apply Dirichlet boundary conditions on the Boundary ids 1 and 2.
Discretization section allows flexibility to change the Finite Element system to be used.
Fuel cell data section is required to set the operating conditions.
Cathode gas diffusion layer is an important section as it determines how the mesh would be treated. We impose that the mesh being imported is a GDL and that the material id 0 is the phase of interest (on which the equation will be solved).
There are options to choose from the different Gas diffusion layer type which would define different values of diffusivity for the simulation from the
OpenFCST database. For our example, we specify a DummyGDL with a diffusion coefficient of 0.22 
.
The next section is Output Variables which is used to specify if the boundary flux is to be computed. By default, the boundary flux is computed on Boundary id 2 for the solution variable which is the molar fraction of oxygen in this case. The final section is Output where the output formats for the solution file can be specified. In addition, to the VTK files, the solver also outputs the solution as .dat file which contains the solution at the nodal points.
Once the two parameter files have been setup, they should be placed into the same folder as the VTK mesh file (${FCST_DIR}/examples/microscale/3D_Diffusion) and the following command can be used to run the simulation:
$ fuel_cell-3d.bin main.prm
If the simulation runs successfully, the solution files called “fuel_cell_solution_DataFile_00001_Cycle_*.vtu” would be written in the same folder. These files can be viewed using Paraview for further postprocessing. The figure below shows a screenshot of the result,
Oxygen molar fraction profile in the mesh after diffusion simulation
Additionally, the oxygen molar flux is also computed and printed to the shell as,
Response oxygen_molar_fraction is : 3.68304e-12
This is an important piece of information as it can be used to compute the effective diffusivity of the media using the following relation:
where 
is the flux output as shown above, 
is the diffusion length along the concentration gradient 
and 
is the cross-section area.
Now, to setup a simulation with Knudsen effects using the local pore radius (as described in the previous section), we have to make some modifications to the main.prm file. In order to include Knudsen effects an additional option called simulator specification = knudsen must be specified. The modified main.prm file is shown below,
###################################################################### # $Id: main.prm 2011-02-09 secanell $ # # # # This file is used to simulate app_diffusion and obtain a # # a concentration profile # # # # # # Copyright (C) 2011 by Marc Secanell # # # ###################################################################### ###################################################################### subsection Simulator set simulator name = diffusion set simulator specification = knudsen set simulator parameter file name = data.prm set nonlinear solver name = None set Analysis type = Analysis end ######################################################################
The data.prm file remains the same as above because computation of the Knudsen effects is done internally and no modifications are required from the user in the data.prm file.
In order to run a simulation with Knudsen, the following command can be executed from the command prompt in the working directory ${FCST_DIR}/examples/microscale/3D_Diffusion
$ fuel_cell-3d.bin main_with_knudsen.prm
Once the simulation runs through the oxygen molar flux at the outlet will be calculated and printed to the shell.
Response oxygen_molar_fraction is : 7.23639e-13
Note
The oxygen flux at the outlet has decreased due to the additional Knudsen resistance introduced in the simulation.
7.2.4. Setting up a diffusion simulation with reaction boundary condition¶
The reaction boundary condition described earlier can be used to simulate oxygen diffusion in a microstructure with reaction at the solid-pore interface. A thin fictitious ionomer film is introduced to account for local losses. Details of the model can be found in [2].
The main.prm file for this case looks as follows,
######################################################################
# $Id: main.prm 2011-02-09 secanell $ #
# #
# This file is used to simulate app_diffusion and obtain a #
# a concentration profile #
# #
# #
# Copyright (C) 2011-15 by ESDLab #
# #
######################################################################
######################################################################
subsection Simulator
set simulator name = diffusion
set simulator specification = reaction #reaction|reaction_and_knudsen|knudsen
set simulator parameter file name = data.prm
set nonlinear solver name = Picard
set Analysis type = Analysis
end
######################################################################
The simulator specification = reaction for this case. This will simulate the reaction BC without Knudsen effects. Additionally, the simulator specification can be changed to reaction_and_knudsen (as shown in the comment) which would simulate the reaction BC and take into account the local Knudsen effects. As described earlier, that this is a non-linear problem therefore, a non-linear solver must be selected. The option nonlinear solver name = Picard specifies the use of Picard iterations.
Several changes need to be made to the data.prm file to setup the additional parameters for the simulation. The updated data.prm file looks as follows,
###################################################################### # $Id: $ # # This file is used to simulate a diffusion application # on a 3D Mesh. # Copyright (C) 2015 by Marc Secanell # ###################################################################### ###################################################################### include ../template/data.prm ###################################################################### ###################################################################### subsection Discretization set Boundary fluxes = true end ###################################################################### ###################################################################### subsection Picard set Under-relaxation = true set Alpha = 6 set Gamma min = 0.3 set Max steps = 10000 set Absolute tolerance = 1e-7 set Relative tolerance = 1e-5 end ###################################################################### ###################################################################### subsection Equations subsection Ficks Transport Equation - oxygen subsection Initial data set oxygen_molar_fraction = 0: 0.1 # where 0 indicates the material_id setup in the grid and 100 is the concentration of solute in mol/cm^3 end subsection Boundary data set species_flux = 0:1 set oxygen_molar_fraction = 1: 0.161, 2: 0.161, 3: 0.161, 4: 0.161, 5: 0.161, 6: 0.161 # where 1 - 6 denote the boundary ids and 0.161 is oxygen molar fraction at these boundary ids end subsection Reaction data set Electronic potential,[V] = 0.2 set Protonic potential,[V] = 0 set Type of kinetics = DoubleTrapKinetics set Reaction name = ORR set Ionomer dissolution constant,[cm/s] = 10 set Ionomer film thickness,[cm] = 10e-7 set A_Pt,s|g = 20 end end end ###################################################################### ###################################################################### subsection Output Variables set Compute boundary responses = true set Output boundary id = 0 set num_output_vars = 1 set Output_var_0 = oxygen_molar_fraction end ######################################################################
The first line in the data.prm includes the template. Therefore, all the parameters presented in the previous section are included and then only the parameters which are required for this example are changed. Since the reaction occurs at the solid-pore boundary and the reaction term is a boundary flux Boundary fluxes parameter in the Discretization section is set to true.
As defined in the main.prm the Picard solver is used to solve the non-linear problem. The section Picard in the parameter file defines some of parameters for under-relaxation and convergence tolerance for the solution variable.
The most important section is the Equations->Ficks Transport Equation - oxygen where the necessary parameters for the reaction boundary condition need to be defined. The Initial data subsection is important for this problem because it is a non-linear problem so a good initial guess is important for the solution. The set oxygen_molar_fraction = 0:0.1 sets the initial oxygen molar fraction in material id 0 which is the pore space for the mesh to 0.1. The Boundary data subsection has two entries. The set oxygen_molar_fraction = 1: 0.161, 2: 0.161, 3: 0.161, 4: 0.161, 5: 0.161, 6: 0.161 sets the oxygen molar fraction boundary condition on the boundary ids 1-6. The set species_flux = 0:1 sets boundary id 0 (which is the boundary id for the solid-pore interface; default value by deal.II) as the boundary where the reaction flux is applied. Reaction data section sets the electronic and protonic potentials for the simulations. The kinetics method is for computation of the current density is selected as DoubleTrapKinetics which is a multi-step model for the oxygen reduction reaction (ORR). Ionomer dissolution constant and Ionomer film thickness define the fictitious ionomer properties at the interface and A_Pt,s|g is the ratio of active area of platinum to the interface area of the solid-void interface.
Output Variables section is modified partially to change the Output boundary id to 0 because we want the flux of oxygen reacting at the solid-pore interface as the output to compute the total current from the simulations.
With these paramters when the simulation is run the output response is,
Response oxygen_molar_fraction is : 1.57756e-14
This translates to a total current of 6.08e-9 A using Faraday’s Law. The figure below shows the oxygen profile screenshot from the simulation,
Oxygen molar fraction profile in the mesh after diffusion with reaction simulation
7.2.5. Additional examples¶
More examples for the microstructures and simulation results reported in Sabharwal et al. [2] are provided in the folder ${FCST_DIR}/examples/microscale/Sabharwal_Analysis_of_FC_CL_2016/. The folder contains two folders namely, Diffusion and Reaction.
Diffusion folder contains the parameter files to simulate the diffusion with and without Knudsen effects and is similar to the example above. The meshes for Stack 1 and Stack 2 from [2] can be downloaded from the OpenFCST website.
Reaction folder contains the parameter files used for the section 3.3.1 for the publication [2]. In order to run the entire polarization curve the variable Equations->Ficks Transport Equation - oxygen->Reaction data->Electronic potential should be varied from 0.9-0.25 V and the output can be converted to the current density using the data provided in [2]. The mesh used for section 3.3.1 titled Mesh1.vtk can be downloaded from the OpenFCST website.
7.3. Electron transport simulation in OpenFCST¶
Similar to the diffusion application electron transport can be simulated in the solid-phase of the microstructures using OpenFCST application AppOhmic. Please check the AppOhmic documentation and examples page for more information on running those simulations.
References
| [1] |
|
| [2] |
|