EECNetworkImpedance.jl

Introduction

The goal is to simulate an impedance measurement of an electrochemical cell of defined structure consisting of YSZ, LSM and pores using electrical elements.

Geometry

The structure is defined by a material matrix which contains integers

• 0 : for a YSZ particle (yellow),
• 1 : for a LSM particle (black),
• 2 : for a pore (white).

Material matrix can be also defined by a raster image using colors specified above in brackes.

In addition, left and right sides of the matrix is grasped by 1 column of LSZ layer as current collectors.

Electrochemistry

Physical simulation is done using standard electrical elements - R elements for ohmic resistance and RC element for a double-layer behavior on interfaces of LSM | YSZ. We have the following parameters:

• R_YSZ : ohmic resistance of YSZ
• CpolYSZ : polarization capacitance on YSZ side
• RpolYSZ : polarization resistance on YSZ side
• R_LSM : ohmic resistance of LSM
• CpolLSM : polarization capacitance on LSM side
• RpolLSM : polarization resistance on LSM side
• R_pore : resistance of pore (is set to very high number as default)

Geometry specific behavior is simulated via interactions between pixels by the manner described in the following picture.

The each impedance Z is specified using the information about starting material M1 and the ending material M2 (see the picture) in the following manner:

• M1 = YSZ
  • M2 = YSZ => Z = R_YSZ/2
  • M2 = LSM => Z = RYSZ/2 + ZRC(RpolYSZ, CpolYSZ) (including the double-layer on YSZ side)
  • M2 = pore => Z = R_YSZ/2
• M1 = LSM
  • M2 = YSZ => Z = RLSM/2 + ZRC(RpolLSM, CpolLSM) (including the double-layer on LSM side)
  • M2 = LSM => Z = R_LSM/2
  • M2 = pore => Z = R_LSM/2
• M1 = pore => Z = R_pore/2

The factor 1/2 is there so that the total resistance through the whole (e. g.) YSZ particle is R_YSZ.

Installation

The package can be then installed via julialang ] add https://github.com/Masicko/EECNetworkImpedance.jl

Usage

Before using the package, you have to execute

julialang using EECNetworkImpedance

Basics

Supposing we have either

• material_matrix = [1 1 1; 0 1 2](with values in {0,1,2}) or
• bitmap image my_image.png (with colors {yellow, black, white}).

In addition, we can specify parameters as a set of pairs. The following are the default parameters:

julialang physical_parameters = [ "R_YSZ" => 100, "R_pol_YSZ" => 0, "C_pol_YSZ" => 0.001, # "R_LSM" => 1, "R_pol_LSM" => 40, "C_pol_LSM" => 0.005, # "R_pore" => 1000000]

Note that "R_pol_YSZ" => 0 which means only one significant arc should appear in Nyquist diagram and therefore two point extrapolation (see below in section Additional options) is viable as resulting to a dramatic speedup. To be precise, more arcs in Nyquist diagram can occure depending on the geometry. For example, even three equally large arcs can be produced using a artificially tailored geometry. However, for an efficient computation on sufficiently large and homogenous geometries, an assumption of one significant arc (and two point extrapolation) is valid.

If less parameters are specified, the others are supposed to be default, i.e.

julialang physical_parameters = ["R_YSZ" => 73]

The core function is

julialang f_list, Z_list = image_to_EIS(material_matrix, physical_parameters)

or using path to image file

julialang f_list, Z_list = image_to_EIS("images/geometry.png", physical_parameters)

or specifying parameters in a function call

julialang f_list, Z_list = image_to_EIS("images/geometry.png", ["R_YSZ" => 73])

or without specifying parameters

julialang f_list, Z_list = image_to_EIS("images/geometry.png")

which returns (by default) frequencies f_list for which impedances Z_list are computed.

Additinal options

Practically useful keyword parameters are

• f_list = [1, 10, 100] : specification of array of frequencies for which EIS simulation will run.
  • = [10.0^n for n in (-3 : 0.5 : 7)] is a good format for a geometric sequence with base 10.0 (it must be floating point number, not integer 10) and exponent spanning from -3 to 7 with a step 0.5.
  • = "two_point_extrapolation" : the simulation is run only for TPE_f_list_in yielding two impedances, R-RC circuit is analytically fitted to the two computed impedances. The output Zlist is computed using this R-RC circuit for frequencies in `TPEflistout`.
  • default value is = "TPE" : which is a shortcut for "twopointextrapolation" with the same meaning
• TPE_f_list_in = [1e-3, 1e6] by default
• TPE_f_list_out = [10.0^n for n in (-3 : 0.5 : 7)] by default
• TPE_warning = true : if true and "TPE" = true, a warning about performing TPE is printed to standard output.
• pyplot = true : if false, no Nyquist plot is plotted
• return_R_RC = false :
  • if = true : the output of function image_to_EIS is a tripple (Rohm, Rpol, C_pol) from R-RC circuit
  • if = false : the output is a tuple (f_list, Z_list)
• export_z_file = "" : decides whether a standard file for Z_view is exported
  • default value is = "", which means do nothing
  • if = "some_file.z" : exports to this file
  • if = "!use_file_name" : this option is valid only when the function image_to_EIS was called with a path of image, e. g. "images/geometry.png" and it means that z_file will have a form "images/geometry.z", i. e. changes only the extension to ".z"
• `savealsoimage = ""
  • if "example.png" : the image will be saved with this name
  • if !asZfile : if export_z_file != "" and than the input image is copied with a name of export_z_file but with the extension ".png"
  • if !input : if image_to_EIS was called with a path of image
• store_R_RC : turns on the evaluation of R_RC element from two points of computed impedance and append the output to a specified file
  • if = "" : means do nothing (which is default)
  • if = "storage.txt" : append a line formated as tab tab tab tab . If the function was called with a matrix (not an path to file), "<matrix_input>" is written instead of .
  • if = "storage.csv" : csv extension works too
• return_specific_impedance = true : returns specific impedance in ohms x cm. Also values of fitted R_RC circuit will be returned in specific units.
• L_el_mat = i_LSM : defines material of the left electrode adjacent to the specified material matrix. Possible options are i_LSM (default) or i_YSZ.
• R_el_mat = i_LSM : the same for the right electrode material.

Advanced keyword parameters are

• complex_type = ComplexF64 : changes the data type in which the impedance calculation is performed
• iterative_solver = "auto" : for small problems (under $15^3$ voxels) is used direct solver and iterative for larger ones.
  • if = false : the system of equations is solved by a direct LU solver (julian \ operator)
  • if = true : the system is solved by iterative solver using Biconjugate gradient stabilized method with using Crout version of incomplete LU decomposition as a preconditioner.
• fill_in_ratio = 12: defines expected fill during the incomplete LU desomposition. It is a ratio of non-zero element count of aproximation of L matrix with respect to non-zero element count of original A matrix of the (sparse) linear system. Bigger fill_in_ratio means more realiable convergence of iterative solver but it requires more RAM space, which is the limiting factor for large systems (e.g. $10^6$ voxels and above).
• tau = "auto": a drop criterion parameter in incomplete LU process. It is a absolute (not relative) treshold and elements of size under this treshold are forgotten. If "auto", an adaptive process is run to estimate a feasible value for tau for which the defined fill_in_ratio is achieved.
• compute_tortuosity = true means that tortuosity of the domain will be computed and stored in the output in the column tor.

Lower level API

A user can directly access the computational core using function material_matrix_to_impedance() which accepts arguments these non-keyword arguments

• material_matrix
• physical_parameters

and keyword arguments

• f_list
• complex_type
• iterative_solver
• verbose : if true, every demanding computation step is measured using a macro @time
• return_only_linsys : if true : the linear system of equations is constructed and returned as a tuple $(A, b)$ - matrix $A$ and right-hand-side $b$. The system is not solved.
• tau
• fill_in_ratio

Real Example

julialang image_to_EIS( [1 0 1; 0 1 2], ["R_YSZ" => 73], # export_z_file="test.z", return_R_RC=true, save_also_image="!asZfile", store_R_RC = "" )

Generate random domain with structure

Automated generating of random structure is essential for statistical testing of system behavior. There are a few helping features using two main parameters.

• porosity in [0, 1] : the ratio of pores over the total points (material points + pores) in the picture.
• LSM_ratio in [0, 1] : probability that the material point will be LSM.

Simple homogenous matrix

The simplest example is a homogenous domain of dimensions = (m, n), where m is a number of rows and n a number of columns. Matrix of this type can be constructed via

julialang homogenous_matrix = generate_matrix(dimensions, porosity, LSM_ratio)

Additional stuctural parameters

More specific structure can be defined using parameters which governs a tendency to group pixels together. Namely:

• pore_cavitance = Nothing- a probability that next pore pixel will spawn next to already existing pore (if possible)
  • if = Nothing, generating of random structure is fast because no "grouping algoritm" is involved.
  • if = 0.0, generating of random structure is slower because "grouping algoritm" is invoked
• LSM_cavitance = 0.0- a probability that next LSM pixel will spawn next to already existing LSM pixel (if possible). This keyword argument is working only if porecavitance keyword parameter is set to a real number (i.e. if `porecavitance = Nothing, then parameterLSM_cavitance` has no effect to generated structure).

julialang pseudo_homogenous_matrix = generate_matrix(dimensions, porosity, LSM_ratio, pore_cavitance = 0.4, LSM_cavitance = 0.2 )

Rotation of domain

Generated domain can be rotated by function

julialang rotate_matrix!(domain, step, axis)

which rotates the input domain (and changes it). There is a version (withoug the exclamation mark !) which does not change the original domain and returns the new rotated domain

julialang new_domain = rotate_matrix(domain, step, axis)

The parameters are - step specifing the number of simple rotations. 1 = 90°, 2 = 180°, 3 = 270°, 4 = 360°, 5 = 450°. And -1 = 90° backwards etc... - axis from the set ["x", "y", "z"] specifies the axis of rotation. Note, the current flows along "y" axis, therefore rotation aroud this axis does not change the output impedance.

Inspecting domain properties

Voxels in domain can be marked by a number, which will specify its material connectivity character. In particular,

• 0 - pore
• 1 - izolated material (not connected to electrodes)
• 2 - material connected to left electrode
• 3 - material connected to right electrode
• 4 - material connected to both electrodes (only this voxels can participate on electric current flow)

This investigation is dome by function

julialang characterize_material_izolation(material_matrix)

which returns matrix labeled as described above. In order to compute ratio of izolated material, there is a function

julialang izolated_material_ratio(material_matrix, wanted_ids = [1,2,3])

which return the ratio of "wantedmaterialpoints / allmaterialpoints". Defaut value is wanted_ids = [1] for strongly izolated voxels.

Structure using multiple submatrices

For more complicated domains composed of several different homogenous subdomains, there is a possibility to construct appropriate matrix. Suppose we want to construct m x n matrix consisting of 2 different submatrices. First, a list of submatrices must be created such that one submatrix is represented by its location (left upper corner and right lower corner) in the resulting matrix and porosity and LSM_ratio.

```julialang # the structure for each submatrix in the list is >> # # [left upper coord, right lower coord, porosity, LSM_ratio]

submatrix_list = [
[(1,1), (10, 5), 0.2, 0.5], # this is the first submatrix [(1,6), (10, 10), 0.0, 1.0] # this is the second submatrix ]

twosubdomainmatrix = generatematrix(submatrixlist) ```

There can be more subdomains in submatrix_list. Dimensions of the resulting two_subdomain_matrix is computed as en rectangular envelope of all locations in submatrix_list. The submatrices can overlap (in this case, the latter has priority in evaluation of matrix), but every pixel must be covered by the submatrices.

Three column domain template

There is a template using the upper structure of defying submatrices, which generates a domain of 3 columns with defined material specification (porosity and LSM_ratio). In addition, there are contacts (of width 1 and optional height) on the left side providing an interesting distribution of electrical current through the system. The right side consists of a continuous one layer of LSM as a connection to conductive electrolyte. The obligate input parameters are

• LSM_ratio1, LSM_ratio2, LSM_ratio3

Optional parameters (with default values) are

• porosity1=0.5, porosity2=0.5, porosity3=0.5
• positions_of_contacts=[15, 50] : starting row for each LSM contact
• height_of_contacts=5
• column_width=5
• height=70

in this default case, the LSM contacts will be between [15, 20] and [50, 55] pixels of the first column while the whole matrix has height of 70 pixels. The submatrix list can be then obtained by function three_column_domain_template and then inserted to generate_matrix.

```julialang LSMratio1 = 0.2 LSMratio2 = 1.0 LSM_ratio3 = 0.5

templatesubmatrixlist = threecolumndomaintemplate(LSMratio1, LSMratio1, LSMratio1, #
columnwidth = 10, porosity1 = 0.0, porosity3 = 0.2, positionsofcontacts=[20, 45], heightof_contacts=4 )

threecolumnmatrix = generatematrix(templatesubmatrix_list) ```

Material matrix visualization

A material matrix can be saved to a file with specified path using matrix_to_file function. For example

julialang matrix_to_file("images/three_column_domain.png", three_column_matrix)

Enlarge image

The images are often very small (50 pixels). So for a better view, there is a utility to make bigger resolution.

julialang enlarge_image(path; resolution=nothing, resize_factor=nothing, print_bool=false)

Parameters are

• path: path for the source image. The processed (enlarged) image is than saved to "path_resized.XXX"
• resolution = nothing you can choose if you want to specify the target resolution as an integer (it assumes squared pictures)
  • resolution = 100 makes a resized picture with dimensions 100x100 pixels
• resize_factor = nothing or you can choose resizing factor, such that
  • resize_factor = 3 means it will make 150x150 if the input image was 50x50
• print_bool = false ... if true, a name of resized picture is printed to terminal

In-depth surface visualization

The function make_shaded_view takes a 3D structure domain and makes a 2D surface view such that it colors visible material pixels according to their depth in the structure. The basic usage is

julialang make_shaded_view(domain, file_name, grad_depth=15, side=1, only_palette = false, LSM_s = (0.3, 0.3, 0.3), LSM_e = (0.0, 0.0, 0.0), YSZ_s = (1.0, 0.8, 0.0), YSZ_e = (0.5, 0.4, 0.0), pall_func = x->x ) Parameters are

• grad_depth = 15
• only_palette = false: if true, the output image is not the 2D surface view but only the palletes of colors
• LSM_s, LSM_e: starting and ending RGB color for gray LSM palette. Similarly, YSZ_s and YSZ_e for yellow YSZ pallete.
• pall_func = x -> x: user can define a function how pallete color should progress from start to end. The x -> x means linear progress. The input parameter x goes from 0 to 1 and the output of the function is expected also in the interval [0, 1]. There is a predefined symmetrical logaritmic function
  • = sym_log(x, a) which has the value 0.5 for x = 0.5 and behaves as a logarithm in interval [0.0, 0.5] and is point-wise symmetric around the midpoint x,y = (0.5, 0.5). The parameter a governs the "logarithmic steepness" in the beginng. a = 1 is almost linear, a = 1000 is very steep
• side = 1: defines, which side od the cube domain is inspected. In particular, the coding is 1 = do not rotate, 2 = 90° along x, 3 = 180° along x, 4 = 270° along x, 5 = 90° along y, 6 = -90° along y

Matrix compilation of subimages

One can generate a lot of images using this package. There is also a tool to make a matrix of images based on row and column indexes:

julialang subimages_composition(folder_path = "data/set_of_images/", row_numbers = collect(0.2 : 0.1 : 0.7), col_numbers = collect(0.2 : 0.05 : 0.4), template = "img_por#(row)_LSM#(col).png", )

where there parameters are

• template = "img_por#(row)_LSM#(col)" where the sequence #(row) will be replaced by a number from the list row_numbers nad the same holds for #(col). Note that it is not $(col) so that it will not colide with the internal julia string interpolation.
• background_color = (1,1,1) RGB definition of background color. White is default.
• del_width = 10 number of pixels separating subimages
• out_filename = "compilation.png" and output file is located in folder_path

Acknowledgement

This work was supported by the German Research Foundation, DFG project no. FU 316/14-1, and by the Czech Science Foundation, GAČR project no. 19-14244J.