英文:
PDE system with variables of different dimensions - Porous electrode theory with FiPy
问题
我理解，不提问，只返回翻译：
# 我试图解决描述多孔电极中稳态电荷守恒的四个偏微分方程系统。方程如下：
1. ∇ ⋅**𝑖**<sub>1</sub> + ∇ ⋅ **𝑖**<sub>2</sub> = 0，
2. **𝑖**<sub>1</sub> = −𝜎 ∇𝜙<sub>1</sub>，
3. **𝑖**<sub>2</sub> = −𝜅 ∇𝜙<sub>2</sub>，
4. ∇ ⋅ **𝑖**<sub>1</sub> = 𝑖<sub>rxn</sub>，或者替代地 ∇ ⋅ **𝑖**<sub>2</sub> = -𝑖<sub>rxn</sub>。
我想使用 FiPy 在 1D 网格上解这些方程作为耦合的偏微分方程组。
# 以下是系统参数、常量和数值参数的代码…
# PDE 变量和系数…
# 边界条件…
# 解方程的函数…
# 运行方程求解…
# 请问是否有其他疑问？
英文:
I'm trying to solve a system of four PDEs describing steady-state charge conservation in a porous electrode. The equations are
-
∇ ⋅**𝑖**<sub>1</sub> + ∇ ⋅ 𝑖<sub>2</sub> = 0,
-
𝑖<sub>1</sub> = −𝜎 ∇𝜙<sub>1</sub>,
-
𝑖<sub>2</sub> = −𝜅 ∇𝜙<sub>2</sub>,
-
∇ ⋅ 𝑖<sub>1</sub> = 𝑖<sub>rxn</sub>, or alternatively ∇ ⋅ 𝑖<sub>2</sub> = -𝑖<sub>rxn</sub>.
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
∇ ⋅ (−𝜎 ∇𝜙<sub>1</sub>) + ∇ ⋅ (−𝜅 ∇𝜙<sub>2</sub>) = 0 and ∇ ⋅ (−𝜎 ∇𝜙<sub>1</sub>) = i<sub>rxn</sub>
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<sub>1</sub> and i<sub>2</sub> 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 𝜙<sub>2</sub> 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.:
-
∇ ⋅ (−𝜎 ∇𝜙<sub>1</sub>) + ∇ ⋅ (−𝜅 ∇𝜙<sub>2</sub>) = 0,
-
∇ ⋅ (−𝜎 ∇𝜙<sub>1</sub>) = 𝑖<sub>rxn</sub>
At x=0
- 𝜙<sub>1</sub> = 0 at x=0,
- ∇𝜙<sub>1</sub> = -I<sub>cell</sub>/𝜎 (i.e. 𝑖<sub>1</sub>=I<sub>cell</sub>)
- ∇𝜙<sub>2</sub> = 0 (i.e. 𝑖<sub>2</sub>=0)
At x=L
- ∇𝜙<sub>1</sub> = 0 (i.e. 𝑖<sub>1</sub>=I<sub>cell</sub>)
- ∇𝜙<sub>2</sub> = -I<sub>cell</sub>/𝜅 (i.e. 𝑖<sub>2</sub>=I<sub>cell</sub>)
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 𝑖<sub>rxn</sub> is defined by the Butler-Volmer equation and is zero if the overpotential is zero. This is however not super important.
- The gradient of 𝜙<sub>2</sub> 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()
´´´
</details>
# 答案1
**得分**: 1
`i1` 和 `i2` 都是一阶的，这是问题的一部分，但它们也由以下公式定义：
```python
i1 = -sigma*phi1.grad # 多孔电极矩阵电流，A/cm2
i2 = -kappa*phi2.grad # 多孔电极电解质电流，A/cm2
这使它们不能作为解变量使用。
即使将 i1 和 i2 保留为一阶的 CellVariables，ConvectionTerm 也不是以这种方式设计的。
> ∇ ⋅ (−𝜎 ∇𝜙1) + ∇ ⋅ (−𝜅 ∇𝜙2) = 0 和 ∇ ⋅ (−𝜎 ∇𝜙1) = irxn
这是 FiPy 的正确方法。我不知道你尝试了什么，但我会写成以下之一：
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)
两者都可以并且都保持到机器精度。
注意:
eq = (chargeConservationEq & matrixPotEq1)（颠倒第一组方程的顺序）对我来说行不通，至少在使用 PETSc 时如此。我理解原因，但我认为我们已经考虑到了这一点。phi2也必须接地：phi2.constrain(0, where=mesh.facesLeft)
英文:
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)
or
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.phi2must 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.
通过集体智慧和协作来改善编程学习和解决问题的方式。致力于成为全球开发者共同参与的知识库，让每个人都能够通过互相帮助和分享经验来进步。
评论