Coupled PDE system - how to define 1st order gradients, convection terms and coefficients

[Robob1715]

Hi,

I am currently learning FiPy and I am trying to implement the PDE system below. The variables I want to solve for are cOH, cO2, phi2 and i2. Currently, I get the following error message:

TypeError: The coefficient must be rank 0 for a rank 0 solution variable.

I am also not sure how to implement the 1st order gradients. For example dcOH/dx in equation 3. Is cOH.grad sufficient? I tried implementing it as ConvectionTerm(coeff =.....) but I always ended up with an error that the coefficient must be a vector.

Thank you very much for your help!

#  %%
from fipy import *
from fipy import CellVariable, Variable, Grid1D, TransientTerm, DiffusionTerm, ImplicitSourceTerm, LinearLUSolver, Viewer
import math
from fipy.tools import numerix
from fipy.solvers.scipy import LinearPCGSolver, LinearLUSolver
import numpy as np

#%%
# Define constants 

DOH = 1e-3  # [cm2/s] # test
DO2 = 1e-3 # [cm2/s]
t0 = 0.78
eps = 0.4
F = 96485.3321
R = 8.314
T = 298
gamma = 1.5
kappa = 12 
cOH_0 = 7.1e-3
a = 3864 # [cm2/cm3]
io2ref = 1e-11
alpha_an = 1.5
alpha_cat = 0.5
cO2 = 1e-7
cO2ref = 1e-7
phi1 = 0 # arbitrary value 
Uref = 0.3027-0.0983
iapp = 0.200 # A/cm2

# %%
# Mesh
L = 1e-3
nx = 1000
msh = Grid1D(dx=L / nx, nx=nx)
x = msh.cellCenters[0]

# initial conditions

cOH = CellVariable(name='OH concentration', mesh=msh, value= 7.0e-3, hasOld= True)
cO2 = CellVariable(name='O2 concentration', mesh=msh, value= 0., hasOld=True)
i2 = CellVariable(name='pore current', mesh=msh, value= 0.,hasOld=True)
phi2 = CellVariable(name='electrolyte potential', mesh=msh, value= 0.,hasOld=True)

# Boundary Conditions
i2.constrain(0, where=msh.facesLeft)
i2.constrain(iapp, where=msh.facesRight)
cOH.constrain(cOH_0, where=msh.facesRight)
cO2.constrain(0, where=msh.facesRight)
phi2.faceGrad.constrain(-iapp/(kappa*eps**gamma), msh.facesRight)

# Source

j2 = a*io2ref*((cOH/cOH_0)**2*numerix.exp(alpha_an*F*(phi1-phi2-Uref )/(R*T))-(cO2/cO2ref)*numerix.exp(-alpha_cat*F*(phi1-phi2-Uref )))
coeff_Ohm = -2*R*T/(cOH*F)*(1-t0+cOH/(2*cOH_0))
DO2_eff = eps**gamma * DO2
DOH_eff = eps**gamma*DOH
kappa_eff = kappa*eps**gamma

# Equations

eq1 = (TransientTerm(var=cOH, coeff=eps) == DiffusionTerm(var=cOH, coeff = DOH_eff) + (t0-1)/F*cOH.grad)

eq2 =(TransientTerm(var=cO2, coeff = eps) == DiffusionTerm(var=cO2, coeff=DO2_eff) + 1/(4*F)*j2)

eq3 = (i2/(kappa_eff) == -phi2.grad + coeff_Ohm*cOH.grad)

eq4 = (i2.grad == j2)

#%%
vi = Viewer((cOH))

dt = 1e-3

solver = LinearLUSolver()
for t in range(100):
    cOH.updateOld()
    cO2.updateOld()
    phi2.updateOld()
    i2.updateOld()

    for sweep in range(5):
        res1 = eq1.sweep(cOH, dt=dt, solver=solver)
        res2 = eq2.sweep(cO2, dt=dt, solver=solver)
        res3 = eq3.sweep(phi2, dt=dt, solver=solver)
        res4 = eq4.sweep(i2, dt=dt, solver=solver)

        print('res1', res1)
        print('res2', res2)
        print('res3', res3)
        print('res3', res4)
    vi.plot()
# %%

[wd15]

So, eq1 has a vector for the source term. That needs to be a scalar. cOH.grad is a vector as written. You can replace it with cOH.grad[0], for example, which doesn't matter for 1D. Think carefully about this if you solve in 2D.

The next problem will be that -alpha_cat*F*(phi1-phi2-Uref) is too large for computing an exponential. You'll need to check your numbers there.

Also, equations 3 and 4 need to be defined with FiPy terms. Currently, they don't have terms, only variables. It means nothing to FiPy and won't work at all. It's normally possible to rework equations like this so they make more sense in the context of the finite volume method (i.e. include transient terms, implicit source terms and diffusion terms). One possible way (I have no idea if this is optimal) would be to run a derivative through equation 3 and use equation 4 to eliminate i_2 from equation 3. You then have a diffusion therm for phi_2 and looks more sensible for FiPy. Probably worth adding in a TransientTerm as well just to stabilize things if you're having trouble converging.

[Robob1715]

Hi wd15,

thank you very much for your help!

I forgot to divide by RT in -alpha_catF(phi1-phi2-Uref). So it is actually exp(-alpha_catF(phi1-phi2-Uref)/(RT)) which makes it a lot smaller.

Could I also define the equations with terms like shown below? I replaced the gradients with convection terms. Running a derivative through eq 3 could be a bit tricky I think.

j2 = a*io2ref*((cOH/cOH_0)**2*numerix.exp(alpha_an*F*(phi1-phi2-Uref)/(R*T))-(cO2/cO2ref)*numerix.exp(-alpha_cat*F*/(R*T))(phi1-phi2-Uref)))
DO2_eff = eps**gamma * DO2
DOH_eff = eps**gamma*DOH
kappa_eff = kappa*eps**gamma
coeff1 = (t0-1)/F
cOHfac = numerix.exp(-0.68818+118.75*cOH**0.5-1030.5*cOH+4004.7*cOH**1.5)/2

main equations

eq1 = TransientTerm(var=cOH, coeff=eps) == DiffusionTerm(var=cOH, coeff = DOH_eff) + coeff1*i2.grad[0]

eq2 = TransientTerm(var=cO2, coeff=eps) == DiffusionTerm(var=cO2, coeff=DO2_eff) + 1/(4*F)*j2

eq3 =  i2/(kappa_eff) == ConvectionTerm(coeff  = [[-1.]], var = phi2) -2*R*T/(cOH*F)*(1-t0+cOHfac)*cOH.grad[0]

eq4 = ConvectionTerm(var=i2, coeff=[[1.]]) == j2

solver = LinearLUSolver()

for t in range(100):
    cOH.updateOld()
    cO2.updateOld()
    phi2.updateOld()
    i2.updateOld()

    for sweep in range(10):
        res1 = eq1.sweep(var=cOH, dt=dt, solver=solver)
        res2 = eq2.sweep(var=cO2, dt=dt, solver=solver)
        res3 = eq3.sweep(var=phi2, dt=dt, solver=solver)
        res4 = eq4.sweep(var=i2, dt=dt, solver=solver)

Now I just get this error for equation 3

'binop' object has no attribute 'sweep'

[wd15]

Could I also define the equations with terms like shown below?

It's not going to work. ConvectionTerms without DiffusionTerms are not easy to solve. Highly recommend doing the math taking the derivative of equation 3. That will give you diffusion terms and coupled diffusion terms and then you might be able to solve as a large coupled equation. You'll have to check the signs on the diffusion terms, but likely they will be stable.

I replaced the gradients with convection terms. Running a derivative through eq 3 could be a bit tricky I think.

Why is it tricky? I'm assuming you're worried about the last term in equation 3, but you don't need to take the derivative of that whole thing. Just apply the derivative and assume it's a diffusion equation for C_OH with a complicated diffusion coefficient. Eventually, you can solve equations as a coupled system. Also, the first term in equation 3 could be implicit in either C_OH or phi_2 depending on which gives better stability.

[wd15]

Maybe close out the Stackoverflow question since it's being dealt with here.

[guyer]

Echoing something Daniel said earlier, I strongly recommend you figure out what your initial equations are in higher dimensions. It can be really misleading to only work with $\partial/\partial x$; it's easy to mix up scalars and vectors and it's very hard to tell the difference between convection and diffusion terms in that notation.

For instance, if I had to guess, your Eq. (1) is probably $\epsilon\frac{\partial c{OH}}{\partial t} = \epsilon^\gamma\nabla\cdot\left(D{OH}\nabla c_{OH}\right) + \frac{t^0 - 1}{F}\nabla\cdot i_2$ so i2.grad[0] is likely to be very wrong, but you want to be sure.

[wd15]

As pointed out, good idea to rewrite equations with vector notation.

[Robob1715]

Thank you for your suggestions, I am learning a lot.

I rewrote the equations. I realized I can also just leave out the previous equation 4. I tried solving it (see code below) but I got the error "Factor exactly singular".

Should I rewrite the source of jO2 in equations 1 & 2 maybe as an implicit source for cOh and cO2, respectively?

In vector notation:

In FiPy:

eps = 0.4
gamma = 1.5
DOH = 1e-3  # [cm2/s] # test
DO2 = 1e-3 # [cm2/s]
DOH_eff = DOH*eps**gamma
DO2_eff = DO2*eps**gamma 
t0 = 0.78
F = 96485.3321
R = 8.314
T = 298
kappa = 64
kappa_eff = kappa*eps**gamma
cOH_0 = 4.1e-3
a = 3864 # [cm2/cm3]
io2ref = 1e-11
alpha_an = 1.5
alpha_cat = 0.5
cO2ref = 1e-7
phi1 = 0 # arbitrary value 
Uref = 0.3027-0.0983
iapp = 0.200 # A/cm2

# Mesh
L = 1e-3
nx = 1000
msh = Grid1D(dx=L / nx, nx=nx)
x = msh.cellCenters[0]

# initial conditions
cOH = CellVariable(name='OH concentration', mesh=msh, value= 7.0e-3, hasOld= True)
cO2 = CellVariable(name='O2 concentration', mesh=msh, value= 0., hasOld=True)
phi2 = CellVariable(name='electrolyte potential', mesh=msh, value= 0.,hasOld=True)

# Boundary Conditions
cOH.constrain(cOH_0, where=msh.facesRight)
cO2.constrain(0, where=msh.facesRight)
phi2.faceGrad.constrain(iapp/(kappa_eff), msh.facesRight)

# Source
eta = phi1-phi2-Uref 
jO2 = a*io2ref*((cOH/cOH_0)**2*numerix.exp(alpha_an*F*eta/(R*T))-(cO2/cO2ref)*numerix.exp(-alpha_cat*F*eta/(R*T)))
coeff1 = (t0-1)/F
coeff2 = 2*R*T*kappa_eff/(cOH*F)*(1-t0+cOH/(2*cOH_0))

# Equations
eq1 = TransientTerm(var=cOH, coeff = eps) == DiffusionTerm(var=cOH, coeff = DOH_eff) + coeff1*jO2

eq2 =TransientTerm(var=cO2, coeff = eps) == DiffusionTerm(var=cO2, coeff=DO2_eff) + 1/(4*F)*jO2

eq3 = DiffusionTerm(coeff = kappa_eff, var = phi2) + DiffusionTerm(coeff = coeff2, var = cOH) + jO2

# Solve
vi = Viewer((cOH))
dt = 1e-8

eqsum = eq1 & eq2  & eq3

solver = LinearLUSolver()
for t in range(100):
    cOH.updateOld()
    cO2.updateOld()
    phi2.updateOld()

    for sweep in range(10):
       eqsum.sweep(dt=dt,solver=solver)

    vi.plot()

Thanks again for your help.

[guyer]

Your equations are virtually not coupled at all (only through eq3). You may get more stable behavior by solving the equations separately, rather than as eqsum.

Depending on the value of eta, jO2 can be incredibly destabilizing.

Physically, Eq. (4) does not make sense to me applied everywhere. This is a description of the redox current and is generally only valid at interfaces between electrode and electrolyte. Eq. (3) on the other hand is a current continuity equation, which should be valid everywhere, I would think with Eq. (4) as a boundary condition.