python - FiPy で無限の計算を修正するには?

Question

私の問題は非常に単純です。解決したい方程式 (コードの最後のセクション) を取得すると、計算は無限になります。メッシュか段差の問題でしょうか。

したがって、私の問題は、2つの類似した3つの偏微分方程式を解くことです。いくつかのパラメーターを変更しましたが、これは常に無限です。

これが私のコードです：

import numpy as np
from fipy import *
import matplotlib.pyplot as plt


##Parameters of the problem##
#############################

L=100*1e-06
nx=int(L/(10*1e-06))
dx=L/nx

n=2.
RR=8.314
T=273.15+20
F=96485.33289

j0=-100*10.    #A/m2
c0= 0.5*1000   #mol/m3
Dc = 5*1e-010  #m2/s
Da= 5*1e-010   #m2/s
uc=Dc/(RR*T)
ua=Da/(RR*T)
zc=n
za=-n

perm=8.854*1e-012*78.54

#############################
#############################


timeStepDuration =  0.05 * dx**2 / (2 * min([Da, Dc]))
steps = 5

mesh = Grid1D(nx=nx, dx=dx)

x = (np.asarray(mesh.cellCenters[0]))

phi0= j0/(-n**2*F**2*c0*(ua+uc))*x

cc= CellVariable(name="cations concentration", mesh=mesh, hasOld=True, value=c0)
ca= CellVariable(name="anions concentration", mesh=mesh, hasOld=True, value=c0)
phi= CellVariable(name="electric potential", mesh=mesh, hasOld=True, value=phi0)


##Boundary conditions########
#############################

# CL cations flux, left
cc.faceGrad.constrain((-1./Dc)*(j0/(n*F)+zc*uc*F*cc.faceValue*phi.faceGrad), mesh.facesLeft)
# CL cations concentration, right
cc.constrain(c0, mesh.facesRight)

# CL anions flux = 0
ca.faceGrad.constrain((-1./Da)*za*ua*F*ca.faceValue*phi.faceGrad, mesh.facesLeft)
# CL anions concentration, right
ca.constrain(c0, mesh.facesRight)

# CL for phi, left
phi.constrain(+RR*T*np.log(cc.faceValue/c0)/(n*F), mesh.facesLeft)
# CL for phi, right
phi.faceGrad.constrain((-j0/(n*F)+Dc*cc.faceGrad)/(n*F*uc*cc.faceValue),mesh.facesLeft)
# phi.faceGrad.constrain(-j0/((n*F)**2*c0*(ua+uc)), mesh.facesRight)

coeff_cc=zc*uc*F
coeff_ca=za*ua*F

#############################
#############################


##Equations##################
#############################

eq_cc = TransientTerm(coeff=1., var=cc) == DiffusionTerm(coeff=Dc,var=cc) + \
    (zc*uc*F*cc.faceValue*phi.faceGrad).divergence

eq_ca = TransientTerm(coeff=1., var=ca) == DiffusionTerm(coeff=Da,var=ca) + \
    (za*ua*F*ca.faceValue*phi.faceGrad).divergence    

eq_phi = DiffusionTerm(coeff=1., var=phi) == (-F/perm)*(zc*cc+za*ca)

coupledEqn = eq_cc & eq_ca & eq_phi

実際、計算は方程式のセクションで終わりがなく、特に最初の 2 つ (cc と ca) では無限です。そして、私が期待しているのは、cc、ca、および phi の進化を示すグラフです。