RFB model by @WenRui-Lv666

rfb.py

"""
Created on Tue Sep 28 19:55:05 2021
@author: lv
In order to make the programming process easier, the same parameters as in the literature are used:
Sadeghi, M. A., Agnaou, M., Kok, M. et al.
Exploring the Impact of Electrode Microstructure on Redox Flow Battery Performance Using a Multiphysics
Pore Network Model [J]. Journal of The Electrochemical Society, 2019, 166(10): A2121–A2130
The only difference from the literature is that the "front" and "back" are used as the boundary of the
electrolyte inflow and outflow.
"""

import openpnm as op
import numpy as np
import scipy as _sp
import matplotlib.pyplot as plt

# Generating network
x_num = 20
y_num = 20
z_num = 10
unit = 1e-6
net = op.network.Cubic(shape=[x_num, y_num, z_num], spacing=unit)
prs = (net['pore.back'] * net['pore.right'] + net['pore.back']
       * net['pore.left'] + net['pore.back'] * net['pore.top'] + net['pore.back'] * net['pore.bottom']
       +net['pore.front'] * net['pore.right'] + net['pore.front'] * net['pore.left'] 
       +net['pore.front'] * net['pore.top'] + net['pore.front'] * net['pore.bottom']
       +net['pore.left'] * net['pore.top'] + net['pore.left'] * net['pore.bottom']
       +net['pore.right'] * net['pore.top'] + net['pore.right'] * net['pore.bottom'])
prs = net.Ps[prs]
thrts = net['throat.surface']
thrts = net.Ts[thrts]
op.topotools.trim(network=net, pores=prs, throats=thrts)
internal = net.pores(['front', 'back', 'left', 'right', 'top', 'bottom'], mode="not")

# Adding geometry
geo = op.geometry.StickAndBall(network=net, pores=net.Ps, throats=net.Ts)
#add half the surface area of the throat to each of the adjacent pores
neighbor_throats = net.find_neighbor_throats(pores=net.Ps, flatten=False)
reaction_area = [] 
for i in net.Ps:
    reaction_area_i = geo['pore.area'][i]
    for j in neighbor_throats[i]:
        reaction_area_i = reaction_area_i + 0.5*geo['throat.area'][j]
    reaction_area.append(reaction_area_i)

# The prebuilt water is used, since most of the electrolyte parameters are the same as water,
# only some parameters need to be modified.
electrolyte = op.phases.Water(network=net)
electrolyte['pore.viscosity'] = 1e-3
electrolyte['throat.viscosity'] = 1e-3
electrolyte['pore.diffusivity'] = 1.15e-9
electrolyte['throat.diffusivity'] = 1.15e-9
#Adding physics
phys = op.physics.GenericPhysics(network=net, phase=electrolyte, geometry=geo)
phys.regenerate_models()

#Performing Stokes flow
flow = op.models.physics.hydraulic_conductance.hagen_poiseuille
phys.add_model(propname='throat.hydraulic_conductance',
               pore_viscosity='pore.viscosity',
               throat_viscosity='throat.viscosity',
               model=flow, regen_mode='normal')
sf = op.algorithms.StokesFlow(network=net, phase=electrolyte)
P_in = 1500 # Pa
P_out = 0  # Pa
sf.set_value_BC(pores=net.pores('front'), values=P_in)
sf.set_value_BC(pores=net.pores('back'), values=P_out)
sf.run()
electrolyte.update(sf.results())

#First iteration
#set relaxation factor w
w = 0.7
F = _sp.constants.physical_constants["Faraday constant"][0]
Vcell = 0.9
#initial guess
V_electrolyte  = []
for i in net.Ps:
    guess = 0.05
    V_electrolyte.append(guess)
Voc = 1.098
res_max = 1

#Continue to iterate
while res_max >0.01:
    op_guess = [Vcell-i-Voc for i in V_electrolyte]
    #Performing advection-diffusion with butler_volmer_conc source term
    diffusive_conductance = op.models.physics.diffusive_conductance.ordinary_diffusion
    phys.add_model(propname='throat.diffusive_conductance', model=diffusive_conductance)
    ad_dif_conductance = op.models.physics.ad_dif_conductance.ad_dif
    phys.add_model(propname='throat.ad_dif_conductance', model=ad_dif_conductance, s_scheme='powerlaw')
    electrolyte['pore.electrolyte_concentration'] = 0
    net['pore.reaction_area'] = reaction_area
    phys['pore.electrolyte_voltage'] = V_electrolyte
    phys['pore.solid_voltage'] = Vcell
    phys['pore.open_circuit_voltage'] = Voc
    BV_params = {
            "z":2,
            "j0": 0.5,
            "c_ref": 1000,
            "alpha_anode": 0.5,
            "alpha_cathode": 0.5
        }
    BV_c = op.models.physics.source_terms.butler_volmer_conc
    phys.add_model(propname='pore.rxn_BV_c',
                            model=BV_c ,
                            X="pore.electrolyte_concentration", **BV_params)
    ad = op.algorithms.AdvectionDiffusion(network=net, phase=electrolyte)
    ad.set_source(propname='pore.rxn_BV_c', pores=internal)
    inlet  = net.pores('front')
    outlet = net.pores('back')
    ad.set_value_BC(pores=inlet, values=900.0)
    ad.set_outflow_BC(pores=outlet)
    ad.run()

    #There are many calculation modes for the reaction rate, and I am not sure which one should be used here.
    Ri_sum_totalnet = ad.rate(pores=net.Ps , mode='group')
    Ri_eachpore = ad.rate(pores=net.Ps , mode='single')
    Ri_throats = ad.rate(throats=net.Ts , mode='single')
    Ri_throats_sum = ad.rate(throats=net.Ts , mode='group')
    Ri_internalpores_sum = ad.rate(pores=internal , mode='group')
    z = 2 #the number of electrons involved in the cathodic reaction
    Rip_sum = Ri_throats_sum*z*F
    #calculating the ohmic resistance of the membrane
    electrode_thickness = z_num * unit #m
    cross_section = (x_num * unit) * (y_num * unit) #m^2
    membrane_electrical_conductance = 8.3
    Rm = electrode_thickness/(cross_section*membrane_electrical_conductance)
    #calculating average voltage loss across the membrane
    del_mem = Rip_sum *Rm
    #overpotential at the membrane interface
    op_membaine_cal = Vcell-del_mem-Voc

    #Performing Ohmic-Conductionn with butler_volmer_voltage source term
    electrolyte['pore.electrolyte_concentration'] = ad['pore.concentration']
    BV_v = op.models.physics.source_terms.butler_volmer_voltage
    phys.add_model(propname='pore.rxn_BV_v',
                            model=BV_v ,
                            X="pore.electrolyte_voltage", **BV_params)
    electrolyte['throat.electrical_conductivity'] = 33.5
    electrolyte['pore.electrical_conductivity'] = 33.5
    electrical_conductance = op.models.physics.electrical_conductance.series_resistors
    phys.add_model(propname='throat.electrical_conductance', model=electrical_conductance)
    oc_boundary = net.pores(["left", "right", 'front', 'back', 'top'], mode="or")
    oc = op.algorithms.OhmicConduction(network=net, phase=electrolyte)
    oc.settings['conductance'] = 'throat.electrical_conductance'
    oc.settings['quantity'] = 'pore.potential'
    oc.set_source(propname='pore.rxn_BV_v', pores=internal)
    oc.set_value_BC(pores=net.pores('bottom'), values=op_membaine_cal)
    # The boundary of the electrolyte voltage at the membrane interface is del_mem, 
    # but I am not sure whether to use del_mem as the boundary value or op_membaine_cal
    # as the boundary value here.
    oc.set_rate_BC(pores=oc_boundary, rates=0)
    oc.run()
    op_cal = oc['pore.potential']
    a = np.array(op_guess)
    b = np.array(op_cal)
    res = abs(a - b) / abs(a)
    res_max = res.max()
    op_next = b*w +a*(1-w)
    V_electrolyte = Vcell - op_next - Voc
print('success')
electrolyte.update(ad.results())
electrolyte.update(oc.results())
fig = op.topotools.plot_coordinates(net, color=sf['pore.pressure'], 
                                    size_by=net['pore.diameter'], 
                                    markersize=50)
fig = op.topotools.plot_coordinates(net, color=ad['pore.concentration'], 
                                    size_by=net['pore.diameter'], 
                                    markersize=50)
fig = op.topotools.plot_coordinates(net, color=oc['pore.potential'], 
                                    size_by=net['pore.diameter'], 
                                    markersize=50)