1

我有一个由 18 个线性（我假设）PDE 组成的系统。但是，它们是形状非常笨拙的 PDE。每次我尝试使用 FiPy 解决耦合 PDE 的系统时，都会出现错误。我已经自己修复了其中的大部分，但我似乎无法克服这个问题。

首先，您将需要稍后将在程序中添加的方程式列表：

import time
from fipy import Variable, FaceVariable, CellVariable, Grid1D, Grid2D, ExplicitDiffusionTerm, TransientTerm, DiffusionTerm, Viewer,PowerLawConvectionTerm, ImplicitSourceTerm
from fipy.tools import numerix
import math
import numpy as np
import sympy as sp
Temperature = 500.0 #Temperature in Celsius

Temp_K = Temperature + 273.15 #Temperature in Kelvin

Ea = [342,7,42,45,34] #Activation energy in order from the list excel/word file
#Frequency factor ko (Initial k)
k_0 =[5.9 * (10**15), 1.3 * (10**13), 1.0 * (10**12), 5.0 * (10**11), 1.2 * (10**13)]
fk1 = [math.exp((-1.0 * x)/(R*Temp_K)) for x in Ea] #Determines k value at given temperature value (exp(-Ea/R*T))
final_k = [fk1*k_0 for fk1,k_0 in zip(fk1,k_0)] #Multiplys by the inital k value to determine the k value at the given temperature (ko * exp(-Ea/R*T))
     final_kcm = [x for x in final_k]
def rhs(eq):
       eq_list = [-1.0*final_kcm[0]*EDC - final_kcm[1]*EDC*R1 - final_kcm[2]*EDC*R2 - final_kcm[3]*EDC*R4 - final_kcm[4]*EDC*R5 - final_kcm[5]*EDC*R6,
       final_kcm[2]*EDC*R2 - final_kcm[8]*EC*R1 + final_kcm[13]*VCM*R2,
       final_kcm[1]*EDC*R1 + final_kcm[6]*R2*R1 + final_kcm[7]*R1*R3 + final_kcm[8]*R1*EC + final_kcm[9]*R1*C11 + final_kcm[10]*R1*C112 +  final_kcm[12]*R1*VCM + 2.0*final_kcm[20]*R1*C2H2,
       2.0*final_kcm[20]*R2*C2H2,
       final_kcm[15]*R5*VCM,

       return eq_list[eq]

这是我用来生成 PDE 系统的程序的其余部分。

EDC = CellVariable(mesh=mesh, hasOld=True, value=10)
EC = CellVariable(mesh=mesh, hasOld=True, value=0)
HCl = CellVariable(mesh=mesh, hasOld=True, value=0)
Coke = CellVariable(mesh=mesh, hasOld=True, value=0)
CP = CellVariable(mesh=mesh, hasOld=True, value=0)

EDC.constrain(10, mesh.facesLeft)
EC.constrain(0., mesh.facesLeft)
HCl.constrain(0., mesh.facesLeft)
Coke.constrain(0., mesh.facesLeft)
CP.constrain(0., mesh.facesLeft)

nsp =18
u_x = [[ [0,]*nsp for n in range(nsp) ]]
for z in range(nsp):
    u_x[0][z][z] = 1.0

eq0 = TransientTerm(var = EDC) == PowerLawConvectionTerm(coeff = u_x, var = EDC) + ImplicitSourceTerm(rhs(0),var = EDC)
eq1 = TransientTerm(var = EC) == PowerLawConvectionTerm(coeff = u_x, var = EC) + ImplicitSourceTerm(rhs(1),var = (EC))
eq2 = TransientTerm(var = HCl) == PowerLawConvectionTerm(coeff = u_x, var = HCl) + ImplicitSourceTerm(rhs(2),var = (HCl))
eq3 = TransientTerm(var = Coke) == PowerLawConvectionTerm(coeff = u_x, var = Coke) + ImplicitSourceTerm(rhs(3),var = (Coke))
eq4 = TransientTerm(var = CP) == PowerLawConvectionTerm(coeff = u_x, var = CP) + ImplicitSourceTerm(rhs(4),var = (CP))


eqn = eq0 & eq1 & eq2 & eq3 & eq4 

if __name__ == '__main__':
     viewer = Viewer(vars = (EDC,EC,HCl,Coke,CP))
     viewer.plot()

for t in range(1):
    EDC.updateOld()
    EC.updateOld()
    HCl.updateOld()
    Coke.updateOld()
    CP.updateOld()

    eqn.solve(dt=1.e-3)

很抱歉代码很长，但我真的无法以其他方式展示它。无论如何，这是它返回的错误：

文件“C:\Users\tjcze\Anaconda3\lib\site-packages\fipy\matrices\scipyMatrix.py”，第 218 行，在 addAt assert(len(id1) == len(id2) == len(vector))

断言错误

我应该怎么做才能让这个系统正常工作？

4

1 回答 1

1

该错误是因为您尝试使用u_x. u_x应该是 rank-1 FaceVariable。它的形状应该是#dims x #faces（或#dims x 1）；它不应该是（#eqns x #eqns）。

设置u_x = [1.]摆脱了AssertionError.

然后，您将收到一系列警告：

UserWarning: sweep() or solve() are likely to produce erroneous results when `var` does not contain floats.

通过初始化所有的CellVariables浮点数来解决这个问题，例如，

EDC = CellVariable(mesh=mesh, hasOld=True, value=10.)

代替

EDC = CellVariable(mesh=mesh, hasOld=True, value=10)

通过这些更改，代码运行。它没有做任何有趣的事情，但这并不奇怪，因为在这个阶段它太复杂了。18 个方程只会混淆事物。我强烈建议您使用 <= 2 个方程来解决此问题。

此时，您的方程根本没有耦合（eq0仅隐式依赖于EDC，eq1仅依赖于EC等）。这没有错，但不是非常有用。eq = eq0 & eq1 & ...语法当然没有意义。EDC 是唯一有机会进化的变量，而且它是恒定的，所以它也不会进化。

将来，请提供实际运行的示例（无论如何，直到出现错误为止）。

于 2019-04-27T00:46:07.430 回答