Как исправить бесконечные вычисления на FiPy? - PullRequest
Новая
       16

Как исправить бесконечные вычисления на FiPy?

0 голосов
/ 30 апреля 2019

Моя проблема довольно проста: когда я получаю уравнения (последний раздел кода), которые я хочу решить, вычисления бесконечны.Может быть, это проблема сетки или шагов?

Так что моя проблема состоит в том, чтобы решить 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.

1 Ответ

0 голосов
/ 01 мая 2019

Этот код входит в бесконечный цикл, потому что ограничения 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)
...