0

我的问题很简单：当我得到我想要解决的方程（代码的最后一部分）时，计算是无止境的。也许是网格或步骤的问题？

所以我的问题是解决3个PDE，有两个相似。我已经更改了一些参数，但这总是无止境的。

这是我的代码：

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

实际上，方程部分的计算是无穷无尽的，尤其是前两个（cc 和 ca）。我期望的是显示 cc、ca 和 phi 演变的图表。

4

1 回答 1

0

此代码进入无限循环，因为 、 和 的约束cc.faceGrad是ca.faceGrad循环phi.faceGrad的。特别是，phi.faceGrad.constrain似乎只是 的算术反转cc.faceGrad.constrain。FiPy 无法解决这个问题。

作为第一步，删除这些约束并将变量与您的方程耦合：

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

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

eq_phi = DiffusionTerm(coeff=1., var=phi) == \
    ImplicitSourceTerm(coeff=(-F/perm)*zc, var=cc) \
    + ImplicitSourceTerm(coeff=(-F/perm)*za, var=ca)
于 2019-05-01T17:21:58.550 回答