Issues with advection diffusion and reactive simulations

[rikvg]

Dear all,

I am currently learning to use OpenPNM for my research on flow batteries and would like to setup a model similar to the one described in:

M. A. Sadeghi, J. Gostick et al., “Exploring the Impact of Electrode Microstructure on Redox Flow Battery Performance Using a Multiphysics Pore Network Model,” J. Electrochem. Soc., vol. 166, no. 10, pp. A2121–A2130, 2019.

For this, I am trying to run advection-diffusion simulations as well as reactive transport simulations. However, I run into an issue when running the code that I do not know how to solve.

First of all, when I try to add the ad_dif_conductance model to a network that has two distinct geometries, the model gives the following error while computing:

File "C:\Users\...\Anaconda3\lib\site-packages\openpnm\models\physics\ad_dif_conductance.py", line 193, in generic_conductance
    Qij = -gh*_sp.diff(P[cn], axis=1).squeeze()

ValueError: operands could not be broadcast together with shapes (1300,) (300,)

The total amount of throats in my system is 1300, however. the geometry that belongs to the physics model that I would like to add the model to only has 300 throats. I believe that I should specify that only the pressure in the pores of geometry type 1 should be used for the calculation of Qij. How would I be able to resolve this error?

Secondly, a question that I have for future implementation: if I want to run a reactive transport simulation, of which I do not know the conversion, should I use the set_rate_BC as output boundary condition by simply setting it to 0 (i.e. is this the same as dC/dx =0)?

Hereby, both scripts (both for the network generation and for the advection-diffusion simulation test (error occurs in the second script in line: phys_big.add_model(propname='throat.ad_dif_conductance', model=mod, s_scheme='powerlaw') ).

Network structure code

import openpnm as op
import openpnm.topotools as tt
import matplotlib.pyplot as plt
import numpy as np

# Clear Workspace
op.Workspace().clear()
plt.close()

# Variables
spacing_distance = 100e-6  # [m] = 100 [um], distance between the big pores
offset_distance = 30e-6  # [m] = 30 [um], distance between big pores and smaller pores

# Starting network for big pores
net = op.network.Cubic(shape=[5, 5, 5], spacing=spacing_distance)

# Save original pores and throat indices and coords
Ps_big = net.Ps
Ts_big = net.Ts

# Make copies of the coords and create new connections
coords = net['pore.coords']
signs = [[-1, -1, 0], [-1, 0, 0], [-1, 1, 0], [0, 1, 0], [1, 1, 0], [1, 0, 0], [1, -1, 0], [0, -1, 0]]  # offset directions

for xsign, ysign, zsign in signs:
    adj = np.zeros_like(coords)
    adj[:, 0] = offset_distance * xsign
    adj[:, 1] = offset_distance * ysign
    adj[:, 2] = offset_distance * zsign
    new_coords = coords + adj
    tt.extend(network=net, pore_coords=new_coords)
    new_Ps = net.Ps[-len(Ps_big):]
    new_Ts = np.vstack((Ps_big, new_Ps)).T
    tt.extend(network=net, throat_conns=new_Ts)

# Check no problems
print(net.check_network_health())

# Save new pores and throat indices
Ps_small = net.Ps[len(Ps_big):]
Ts_small = net.Ts[len(Ts_big):]  # gives all throats from big to small pores

# Apply geometries
geo_big = op.geometry.GenericGeometry(network=net, pores=Ps_big, throats=Ts_big)
geo_small = op.geometry.GenericGeometry(network=net, pores=Ps_small, throats=Ts_small)

# Assign pore and throat properties

big_size = 20e-6
small_size = 2e-6

geo_big['pore.diameter'] = big_size
geo_small['pore.diameter'] = small_size

geo_big['throat.diameter'] = big_size
geo_small['throat.diameter'] = small_size

geo_big.add_model(propname='pore.volume', model=op.models.geometry.pore_volume.sphere)
geo_small.add_model(propname='pore.volume', model=op.models.geometry.pore_volume.sphere)

# geo_big.add_model(propname='throat.diameter', model=op.models.geometry.throat_size.from_neighbor_pores)
# geo_small.add_model(propname='throat.diameter', model=op.models.geometry.throat_size.from_neighbor_pores, mode='min') # Only necessary if smaller pores are interconnected

geo_big.add_model(propname='throat.endpoints', model=op.models.geometry.throat_endpoints.spherical_pores)
geo_small.add_model(propname='throat.endpoints', model=op.models.geometry.throat_endpoints.spherical_pores)

geo_big.add_model(propname='throat.area', model=op.models.geometry.throat_area.cylinder)
geo_small.add_model(propname='throat.area', model=op.models.geometry.throat_area.cylinder)

geo_big.add_model(propname='pore.area', model=op.models.geometry.pore_area.sphere)
geo_small.add_model(propname='pore.area', model=op.models.geometry.pore_area.sphere)

geo_big.add_model(propname='pore.surface_area', model=op.models.geometry.pore_surface_area.sphere)
geo_small.add_model(propname='pore.surface_area', model=op.models.geometry.pore_surface_area.sphere)

geo_big.add_model(propname='throat.surface_area', model=op.models.geometry.throat_surface_area.cylinder)
geo_small.add_model(propname='throat.surface_area', model=op.models.geometry.throat_surface_area.cylinder)

geo_big.add_model(propname='throat.conduit_lengths', model=op.models.geometry.throat_length.conduit_lengths)
geo_small.add_model(propname='throat.conduit_lengths', model=op.models.geometry.throat_length.conduit_lengths)

# Plot
fig = tt.plot_connections(network=net, throats=Ts_big, c='b')
fig = tt.plot_connections(network=net, throats=Ts_small, fig=fig, c='y')
fig = tt.plot_coordinates(network=net, pores=Ps_big, c='r', s=20, fig=fig)
fig = tt.plot_coordinates(network=net, pores=Ps_small, c='g', s=2, fig=fig)

plt.show()

# Save generated structure
op.io.OpenpnmIO.save_project(project=net.project, filename='./output/bimodal_geometry_%d_surrounding_small_pores' %len(signs))
op.io.VTK.save(network=net, filename='./output/bimodal_geometry_%d_surrounding_small_pores' %len(signs)) 

Network Simulation code:

import openpnm as op
import openpnm.topotools as tt
import matplotlib.pyplot as plt
import numpy as np

# Clear Workspace
op.Workspace().clear()

file = './output/bimodal_geometry_8_surrounding_small_pores'

# Load in project:
project = op.io.OpenpnmIO.load_project(filename=file)
net = project.network
geo_big = project['geo_01']
geo_small = project['geo_02']

# Add phases:
water = op.phases.Water(network=net)

# Create physics object for each geometry:
phys_big = op.physics.Standard(network=net, phase=water, geometry=geo_big)
phys_small = op.physics.Standard(network=net, phase=water, geometry=geo_small)

# Run StokesFlow to obtain the pressure field and velocity field in the pore network
P1 = 202650  # Pa
P2 = 101325  # Pa
delta_P = P1 - P2

sf = op.algorithms.StokesFlow(network=net, phase=water)
sf.set_value_BC(values=P1, pores=net.pores('left'))
sf.set_value_BC(values=P2, pores=net.pores('right'))
sf.run()

water.update(sf.results())

mod = op.models.physics.ad_dif_conductance.ad_dif
phys_big.add_model(propname='throat.ad_dif_conductance', model=mod, s_scheme='powerlaw')
phys_small.add_model(propname='throat.ad_dif_conductance', model=mod, s_scheme='powerlaw')

ad = op.algorithms.AdvectionDiffusion(network=net, phase=water)
# ad.settings.update({'s_scheme': 'exact'})

ad.set_value_BC(pores=net.pores('left'), values=100.0)
ad.set_value_BC(pores=net.pores('right'), values=0.0)
ad.run()

#example = op.algorithms.ReactiveTransport

op.io.OpenpnmIO.save_project(project=net.project, filename=file + '_simulated')
op.io.VTK.save(network=net, phases=water, filename=file + '_simulated')

Thank you, Rik

[jgostick]

It might take a while to work out an answer to your question. We're busy this week getting ready for a conference :-)

[rikvg]

Dear Prof. Gostick,

Thank you for the quick response. I will try to work on a solution myself in the meantime, and I am interested in your response when you have the time.

Good luck with your preparations for the conference!

[rikvg]

Dear Prof. Gostick, @magnaou and team,

I have worked further on the topic and decided to first focus on getting the model working for a well-defined cubic network :). I would like to use the ChargeConservationNernstPlanck algorithm to simulate steady-state reactive transport together with ion transport for the vII/vIII redox couple, but I still have some questions about its implementation that I was hoping you might be willing to answer.

The set-up of my code is as following:

• Network generation & phase, geometry and physics instantiation
• Run the Stokes flow algorithm to solve the pressure field
• Run the advection-diffusion model together with a reactive term (Butler Volmer-equation) with initial concentrations to get an idea of the conversion for a given overpotential
• Calculate the estimated equilibrium potential and subsequently the potential field by: vanadiumII['pore.potential_field'] = eta - potential_solid - equilibrium_potential

After this I would like to use the coupled solver as given in the ChargeConservationNernstPlanck algorithm in order to solve the steady state charge transfer and ion transfer.

However, when I call the algorithm, I get the following error:

cc_vII = op.algorithms.ChargeConservationNernstPlanck(network=net, phase=vanadiumII)
cc_vII.settings['ions'] = vanadiumII
cc_vII.settings['potential_field'] = vanadiumII['pore.potential_field']
cc_vII.run()

File "C:\Users\...\Anaconda3\lib\site-packages\openpnm\algorithms\ChargeConservationNernstPlanck.py", line 53, in run
    p_alg[p_alg.settings['quantity']]
AttributeError: 'numpy.ndarray' object has no attribute 'settings'

How would I go about correctly calling this algorithm? Furthermore, the algorithm calls _run_reactive for both the potential field and the concentration of the ions. What equations are exactly solved for both cases? Where would I be able to alter these equations so that I can play with the transport and reactive terms that are taken into account?

Thank you in advance!

Kind regards, Rik