PDE system with variables of different dimensions - Porous electrode theory with FiPy

I'm trying to solve a system of four PDEs describing steady-state charge conservation in a porous electrode. The equations are

∇ ⋅𝑖₁ + ∇ ⋅ 𝑖₂ = 0,
𝑖₁ = −𝜎 ∇𝜙₁,
𝑖₂ = −𝜅 ∇𝜙₂,
∇ ⋅ 𝑖₁ = 𝑖_rxn, or alternatively ∇ ⋅ 𝑖₂ = -𝑖_rxn.

I want to use FiPy to solve these equations as a set of coupled PDEs on a 1D grid.

import fipy.tools.numerix as nrx
from fipy import Grid1D, CellVariable, PowerLawConvectionTerm, DiffusionTerm, Matplotlib1DViewer


## SYSTEM PARAMETERS AND CONSTANTS

cd = 4        # operating current density, A/cm2
T = 333       # operating temperature, K
L = 10*1e-4   # electrode thickness, cm
iEx = 1e-2    # exchange current density, A/cm2
Er = 0        # reaction activation energy, V
nr = 4        # number of electrons in reaction
F = 96485     # faraday's constant, C/mol
R = 8.3145    # molar gas constant, J/Kmol


## NUMERICAL PARAMETERS

nx = 100                          # number of mesh points
dx = L/nx                         # distance between mesh points
mesh = Grid1D(nx = nx, dx = dx)   # 1D grid of length L
tol = 1e-3                        # tolerance for residual
itmax = 100                       # max no. of sweeps to break convergence loop


## PDE VARIABLES AND COEFFICIENTS

irxn = CellVariable(mesh = mesh, value = iEx, name = r'$i_\mathrm{rxn}$')  # reaction current density (i_rxn), A/cm2

sigma = CellVariable(mesh = mesh, value = 50, name = r'$\sigma$')       # electron conductivity, S/cm
kappa = CellVariable(mesh = mesh, value = 100e-3, name = r'$\kappa$')   # proton conductivity, S/cm

phi1 = CellVariable(mesh = mesh, name = r'$\Phi_1$')   # porous electrode matrix potential, V
phi2 = CellVariable(mesh = mesh, name = r'$\Phi_2$')   # porous electrode electrolyte potential, V
i1 = -sigma*phi1.grad                                  # porous electrode matrix current, A/cm2
i2 = -kappa*phi2.grad                                  # porous electrode electrolyte current, A/cm2
i1.name = r'$i_1$'
i2.name = r'$i_2$'


## BOUNDARY CONDITIONS

phi1.constrain(0, where=mesh.facesLeft)       # matrix potential at x=0
i1.constrain(0, where=mesh.facesRight)        # matrix current at x=L
i2.constrain(0, where=mesh.facesLeft)         # electrolyte current at x=0
i2.constrain(cd, where=mesh.facesRight)       # electrolyte current at x=L


## SOLVE EQUATIONS

def solveEquations():

  # equations

  # charge conservation ∇ ⋅ 𝑖_1 + ∇ ⋅ 𝑖_2 = 0 (as ∇ ⋅ [−𝜎 ∇𝜙_1] + ∇ ⋅ [−𝜅 ∇𝜙_2] = 0)
  #chargeConservationEq = (DiffusionTerm(coeff = -sigma, var = phi1) + DiffusionTerm(coeff = -kappa, var = phi2) == 0) - this is commented out for testing but gives a similar error
  chargeConservationEq = (PowerLawConvectionTerm(coeff = (1,), var = i1) + PowerLawConvectionTerm(coeff = (1,), var = i2) == 0)
  

  # matrix potential ∇ ⋅ 𝑖_1 = 𝑖_rxn
  matrixPotEq = (PowerLawConvectionTerm((1,), var = i1) == irxn)
        
  # electrolyte potential ∇ ⋅ 𝑖_2 = -𝑖_rxn
  electrolytePotEq = (PowerLawConvectionTerm((1,), var = i2) == -irxn)

  # matrix current 𝑖_1 = −𝜎 ∇𝜙_1 (divergence taken on both sides)
  matrixCurrentEq = (PowerLawConvectionTerm((1,), var = i1) == DiffusionTerm(coeff = -sigma, var = phi1))

  # electrolyte current 𝑖_2 = −𝜅 ∇𝜙_2 (divergence taken on both sides)
  electrolyteCurrentEq = (PowerLawConvectionTerm((1,), var = i2) == DiffusionTerm(coeff = -kappa, var = phi2))

  eq = ( chargeConservationEq
    & matrixPotEq
    #& electrolytePotEq
    & matrixCurrentEq
    & electrolyteCurrentEq) 

  res = eq.sweep()
  it = 0

  while ((res > tol) and (it < itmax)):

    # increment
    it += 1

    # solve equations
    res = eq.sweep()


solveEquations()

However when I try to run I get the error message: ValueError: all the input arrays must have same number of dimensions, but the array at index 0 has 2 dimension(s) and the array at index 2 has 1 dimension(s). The error occurs when eq.sweep() is called.

I understand that sigma, kappa, irxn, phi1 and phi2 have shape (100,) while i1 and i2 have shape (1,100). Is this why I get this error? If so, how do I write the equations in a way that can be solved?

I have also tried to combine the equations to get

∇ ⋅ (−𝜎 ∇𝜙₁) + ∇ ⋅ (−𝜅 ∇𝜙₂) = 0 and ∇ ⋅ (−𝜎 ∇𝜙₁) = i_rxn

in order to simplify the problem to only two variables. In this case I did not get this error, but I also did not get a stable solution and for some reason charge was not conserved even though the first equation should ensure this - when plotting the currents i₁ and i₂ their sum was not constant as it should be. The boundary conditions were also less stable as there were boundary conditions on the gradients of 𝜙₂ on either side, following from eq. 3.

Any help is appreciated!

EDIT: a new version of the code implementing the suggested changes, with more math details. This code runs without the error originally mentioned, but some new questions/issues have appeared.

A summary of the governing equations boundary conditions that should be implemented.:

∇ ⋅ (−𝜎 ∇𝜙₁) + ∇ ⋅ (−𝜅 ∇𝜙₂) = 0,
∇ ⋅ (−𝜎 ∇𝜙₁) = 𝑖_rxn

At x=0

𝜙₁ = 0 at x=0,
∇𝜙₁ = -I_cell/𝜎 (i.e. 𝑖₁=I_cell)
∇𝜙₂ = 0 (i.e. 𝑖₂=0)

At x=L

∇𝜙₁ = 0 (i.e. 𝑖₁=I_cell)
∇𝜙₂ = -I_cell/𝜅 (i.e. 𝑖₂=I_cell)

The results here are closer to what I would expect, but something is still not right

No matter how low or high I set the exchange current density, nothing seems to happen to the currents i1 and i2. This should not be the case, as a lower magnitude of the exchange current density should mean less generated current at a given overpotential.
Constraining phi2 to be zero at the x=0 boundary implies that there are no reactions there, as the reaction current 𝑖_rxn is defined by the Butler-Volmer equation and is zero if the overpotential is zero. This is however not super important.
The gradient of 𝜙₂ does not appear to be zero at x=0 as specified?


    from fipy import Grid1D, CellVariable, DiffusionTerm, Matplotlib1DViewer
    
    ## SYSTEM PARAMETERS AND CONSTANTS
    
    cd = 4        # operating current density, A/cm2
    T = 333       # operating temperature, K
    L = 10*1e-4   # electrode thickness, cm
    iEx = 1e-2    # exchange current density, A/cm2
    Er = 0        # reaction activation energy, V
    nr = 4        # number of electrons in reaction
    F = 96485     # faraday's constant, C/mol
    R = 8.3145    # molar gas constant, J/Kmol
    beta = 0.5    # symmetry factor for reaction
    
    
    ## NUMERICAL PARAMETERS
    
    nx = 100                          # number of mesh points
    dx = L/nx                         # distance between mesh points
    mesh = Grid1D(nx = nx, dx = dx)   # 1D grid of length L
    tol = 1e-3                        # tolerance for residual
    itmax = 100                       # max no. of sweeps to break convergence loop
    
    
    ## PDE VARIABLES AND COEFFICIENTS
    
    irxn = CellVariable(mesh = mesh, value = iEx, name = r'$i_\mathrm{rxn}$')  # reaction current density (i_rxn), A/cm2
    
    sigma = CellVariable(mesh = mesh, value = 50, name = r'$\sigma$')       # electron conductivity, S/cm
    kappa = CellVariable(mesh = mesh, value = 100e-3, name = r'$\kappa$')   # proton conductivity, S/cm
    
    phi1 = CellVariable(mesh = mesh, name = r'$\Phi_1$')   # porous electrode matrix potential, V
    phi2 = CellVariable(mesh = mesh, name = r'$\Phi_2$')   # porous electrode electrolyte potential, V
    
    
    ## BOUNDARY CONDITIONS
    
    phi1.constrain(0, where=mesh.facesLeft)                                       # matrix potential grounded at x=0
    phi2.constrain(0, where=mesh.facesLeft)                                       # electrolyte potential grounded at x=0
    phi1.faceGrad.constrain(-cd/sigma.faceValue[0], where=mesh.facesLeft)         # matrix current at x=0
    phi2.faceGrad.constrain([0], where=mesh.facesLeft)                            # electrolyte current at x=0
    phi1.faceGrad.constrain([0], where=mesh.facesRight)                           # matrix current at x=L
    phi2.faceGrad.constrain([-cd/kappa.faceValue[-1]], where=mesh.facesRight)     # electrolyte current at x=L
    
    
    ## SUPPORTING FUNCTIONS
    
    def rcd(phiM, phiE):
    
      # calculate reaction current density
      # inputs: phiM (potential in solid electron-conducting phase) & phiE (potential in proton-conducting electrolyte)
    
      # overpotential
      overPotential = phiM-phiE-Er
    
      # Butler-Volmer equation for reaction current density
      i = iEx*(nrx.exp((beta*nr*F)/(R*T)*overPotential)-nrx.exp(-((1-beta)*nr*F)/(R*T)*overPotential))
    
      return i
    
    
    ## SOLVE EQUATIONS
    
    def solveEquations():
    
      # equations
    
      # charge conservation ∇ ⋅ 𝑖_1 + ∇ ⋅ 𝑖_2 = 0 (as ∇ ⋅ [−𝜎 ∇𝜙_1] + ∇ ⋅ [−𝜅 ∇𝜙_2] = 0)
      chargeConservationEq = (DiffusionTerm(coeff = -sigma, var = phi1) + DiffusionTerm(coeff = -kappa, var = phi2) == 0)
      
      # matrix potential ∇ ⋅ 𝑖_1 = 𝑖_rxn
      matrixPotEq = (DiffusionTerm(coeff = -sigma, var = phi1) == irxn)
            
      eq = (matrixPotEq & chargeConservationEq)
    
      res = eq.sweep()
      it = 0
    
      while ((res > tol) and (it < itmax)):
    
        # update value of irxn based on overpotential from previous iteration
        irxn.value = rcd(phi1.value, phi2.value)
    
        # increment
        it += 1
    
        # solve equations
        res = eq.sweep()
    
      print(it,'iterations')
      print('residual',res)
    
    
    def plotPotentials():
    
      viewer = Matplotlib1DViewer(vars = (phi1, phi2), title = r'Potential [V]', legend = 'best')
      viewer.plot()
      input('Press return to continue')
    
    def plotCurrents():
    
      i1 = -sigma.faceValue*phi1.faceGrad                                  # porous electrode matrix current, A/cm2
      i2 = -kappa.faceValue*phi2.faceGrad                                  # porous electrode electrolyte current, A/cm2
      i1.name = r'$i_1$'
      i2.name = r'$i_2$'
    
      x = mesh.faceCenters.value[0]
      plt.figure()
      plt.plot(x,i1.value[0],label=r'$i_1$')
      plt.plot(x,i2.value[0],label=r'$i_2$')
      plt.legend(loc = 'best')
      plt.show()
    
    
    solveEquations()
    
    ## UNCOMMENT TO PLOT RESULTS
    
    #plotPotentials()
    #plotCurrents()

´´´

Solution

i1 and i2 are both rank-1, which is part of the problem, but they are also defined by formulae:

i1 = -sigma*phi1.grad                                  # porous electrode matrix current, A/cm2
i2 = -kappa*phi2.grad                                  # porous electrode electrolyte current, A/cm2

which makes them unusable as solution variables.

Even if i1 and i2 were left as rank-1 CellVariables, ConvectionTerm is not designed to work this way.

∇ ⋅ (−𝜎 ∇𝜙1) + ∇ ⋅ (−𝜅 ∇𝜙2) = 0 and ∇ ⋅ (−𝜎 ∇𝜙1) = irxn

This is the correct approach for FiPy. I don't know what you tried, but I would write either:

chargeConservationEq = (DiffusionTerm(coeff = -sigma, var = phi1) + DiffusionTerm(coeff = -kappa, var = phi2) == 0) 
matrixPotEq1 = (DiffusionTerm(coeff=-sigma, var = phi1) == irxn)
eq = (matrixPotEq1 & chargeConservationEq)

matrixPotEq1 = (DiffusionTerm(coeff=-sigma, var = phi1) == irxn)
matrixPotEq2 = (DiffusionTerm(coeff=-kappa, var = phi2) == -irxn)
eq = (matrixPotEq1 & matrixPotEq2)

Either works and either conserves charge to machine precision.

Notes:

eq = (chargeConservationEq & matrixPotEq1) (reversing the order of the first set of equations) does not work for me, at least with PETSc. I understand why, but I thought we had accounted for this.
phi2 must also be grounded: phi2.constrain(0, where=mesh.facesLeft)

Edit: Revised for revised final question

You are trying to apply six boundary conditions for two 2nd-order PDEs. This is over-constrained. Moreover, trying to set the value and the gradient at the same boundary amounts to a shooting method, which FiPy doesn't readily deal with.

If I change the boundary conditions to:

phi1.constrain(0, where=mesh.facesRight)                                       # matrix potential grounded at x=L
phi2.constrain(0, where=mesh.facesLeft)                                       # electrolyte potential grounded at x=0
phi1.faceGrad.constrain([-cd/sigma.faceValue], where=mesh.facesLeft)         # matrix current at x=0
phi2.faceGrad.constrain([-cd/kappa.faceValue], where=mesh.facesRight)     # electrolyte current at x=L

then I get

which makes some sort of sense.