How can I use PyBaMM to solve a Poisson equation specifying both boundary conditions as Neumann conditions?

Hy everyone. I’m developing a transient 1D model for the membrane of a vanadium redox flow battery. The membrane is treated as a single spatial domain that is discretised according to the finite volume method:


# Spatial variable

x = pybamm.SpatialVariable("x", domain=["membrane"], coord_sys="cartesian")

# Geometry

geometry = {"membrane": {x: {"min": pybamm.Scalar(0), "max": dm}}}

# Create mesh and discretise

submesh_types = {"membrane": pybamm.Uniform1DSubMesh}

var_pts = {x: 100}

mesh = pybamm.Mesh(geometry, submesh_types, var_pts)

spatial_methods = {"membrane": pybamm.FiniteVolume()}

disc = pybamm.Discretisation(mesh, spatial_methods)

The objective is to model the transport of chemical species across the membrane along the longitudinal (x-coordinate) direction. Right now the model is simplified and accounts for the conservation of chemical species and the conservation of electrical charge, which I am modeling with the following equations:


# Variables

c2 = pybamm.Variable("V2 concentration", domain="membrane")

c3 = pybamm.Variable("V3 concentration", domain="membrane")

c4 = pybamm.Variable("V4 concentration", domain="membrane")

c5 = pybamm.Variable("V5 concentration", domain="membrane")

phi_l = pybamm.Variable("Ionic potential", domain="membrane")

# Nernst-Planck equation (considering only diffusion)

N2 = - D2 * pybamm.grad(c2)

N3 = - D3 * pybamm.grad(c3)

N4 = - D4 * pybamm.grad(c4)

N5 = - D5 * pybamm.grad(c5)

# Species conservation

dc2dt = - pybamm.div(N2)

dc3dt = - pybamm.div(N3)

dc4dt = - pybamm.div(N4)

dc5dt = - pybamm.div(N5)

model.rhs = {c2: dc2dt, c3: dc3dt, c4: dc4dt, c5: dc5dt}

# Charge conservation (Simplified: Ohm's Law)

i_l = - sigma_m * pybamm.grad(phi_l)

charge_conservation = pybamm.div(i_l)

model.algebraic = {phi_l: charge_conservation}

I am trying to solve the Poisson equation for the ionic potential by providing two Neumann boundary conditions that take into account the continuity of the current at the electrode/membrane interface:


phi_l: {

"left": (- J_e / sigma_m, "Neumann"),

"right": (J_e / sigma_m, "Neumann")

}

Note that J_e and sigma_m are known values, and the boundary conditions for the concentrations (c2, c3, c4, c5) are all Dirichlet boundary conditions. If I try to solve this problem I get an error: pybamm.expression_tree.exceptions.SolverError: Could not find consistent states: Could not find acceptable solution: solver returned NaNs

This is reasonable because at the boundaries I’m not fixing the value of phi_l, but its derivative, and so there is a floating potential issue. To solve this “floating” problem I would like to force the potential to assume a specific value at one node. To do this, my idea would be to access the finite difference matrix and modify one row by setting all entries to zero, except the entry related to the node where I want to force the potential to have a specific known value. From what I understand in the source code, in the FiniteVolume() class the gradient() method applies the gradient operator to the discretized symbol taking into account the two Neumann conditions on the boundaries, and subsequently the divergence() method applies the divergence operator on the gradient, basically multiplying the divergence matrix by the result of the gradient() method. Is it possible to directly access the finite difference matrix, i.e. the product between the divergence matrix and the gradient matrix, and modify it by hand? Otherwise, how could I do it?

Do you think it is possible to handle this kind of problem in PyBaMM? Have any of you ever worked on a similar problem with two Neumann boundary conditions?