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()
是一种方法,而不是 属性(它需要括号)。
我希望表示 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()
是一种方法,而不是 属性(它需要括号)。