2D 网状网络 fipy 的质量平衡不正确

incorrect mass balance with 2D mesh network fipy

我希望表示 2D 网络中的扩散(扩散系数取决于 phi 的值)和特定单元格中的设定 phi 输入速率(因此不是脸上的 BC)。这似乎是一个非常简单的场景,但是,我一定是做错了什么,因为在计算这个例子时我得到了非常奇怪的结果:

from fipy import *
from fipy.meshes.nonUniformGrid1D import NonUniformGrid1D as Grid1D
from fipy.meshes.nonUniformGrid2D import NonUniformGrid2D as Grid2D
from fipy.tools import numerix as np

ny = 6 
nx = 1 
dy = [0.2, 0.2, 0.2, 0.2, 0.2, 0.2] 
dx = 0.2 
translate = [[0.0], [-1.2]]
meshes = Grid2D(ny = ny,nx = nx, dy = dy, dx = dx) + translate


nx = 2 
ny = 1 
dx = [0.05, 0.05]
dy = 0.2 
translate = [[0.2], [-0.2]] 
mesh = Grid2D(ny = ny,nx = nx, dy = dy, dx = dx) + translate

meshes = meshes + mesh

nx = 2 
ny = 1 
dx = [0.05 ,0.05]
dy = 0.2 
translate = [[0.2], [-0.6]] 
mesh = Grid2D(ny = ny,nx = nx, dy = dy, dx = dx) + translate

meshes = meshes + mesh

nx = 2 
ny = 1 
dx = [0.05, 0.05]
dy = 0.2 
translate = [[0.2], [-1.0]] 
mesh = Grid2D(ny = ny,nx = nx, dy = dy, dx = dx) + translate

meshes = meshes + mesh

ny = 3 
nx = 1 
dy = [0.2 ,0.2, 0.2]
dx = 0.2 
mesh = Grid2D(ny = ny,nx = nx, dy = dy, dx = dx) 

meshes = meshes + mesh


ny = 1 
nx = 1 
dx = [1.]
dy = 0.2 
translate = [[0.2], [0.2]]
mesh = Grid2D(ny = ny,nx = nx, dy = dy, dx = dx) + translate

meshes = meshes + mesh
vol = meshes.cellVolumes
phi = CellVariable(mesh= meshes, value = 0., hasOld = True)
Source = CellVariable(mesh= meshes, value = 0.)
Source.constrain(10., where = abs(vol - 0.2) < 0.001) 
#addition of Source * dt to meshes at each step on the selected cell
dt = 1e-6#0.9*( min(meshes.cellVolumes) ** 2) / (2 *max(Source.faceValue ))

cumulatedIn = 0. # cumulated input
eq = (TransientTerm(var = phi) == DiffusionTerm(var = phi,coeff= phi.faceValue)+ Source )
steps = 3
for _ in range(steps):
   print('steps ', _, dt)
   res= 1e5
   loop =0 
   while (res > 1e-3* max(abs(Source.value)) or res > 1e-3* max(abs(phi.value)) )and loop < 1000:
           res =  eq.sweep(dt= dt)
           loop += 1
           print('res: ',res)
           print('maxRes: ',1e-3* max(abs(Source.value)) , 1e-3* max(abs(phi.value)) )
   cumulatedIn += sum(Source * dt * meshes.cellVolumes)
   phi.updateOld
   print('mass balance ', sum(phi * meshes.cellVolumes) - cumulatedIn,'\nphi content in mesh ', sum(phi * meshes.cellVolumes))
   print('phi matrix\n', phi)

phi的含量没有增加而是保持不变。 当我改为使用 solve 时,该值确实增加了。但是,因为扩散系数取决于 phi,所以我宁愿使用扫描。 有人可以解释问题的根源吗?

2) 此外,扩散项应该表示浓度驱动的平流:J = [constante coeff to be added later] * nabla_dot_(phi * nabla(phi))。 或者用 fipy 显式方程写成: (constante_coeff * phi * phi.faceGrad).发散。 我当前的扩散术语是表示这个的正确方法吗?

感谢您的帮助!

.updateOld() 是一种方法,而不是 属性(它需要括号)。