fipy - 一般边界条件

Question

我很感激有关设置一般边界条件的帮助，-grad(y) + g(y) = 0其中g是一些未知函数y。这是一个我无法开始工作的简单一维示例：

N=3
h=1./(float(N)-1.)

mesh = Grid1D(nx=N, dx=h)

c=CellVariable(mesh=mesh,value=0.5)

## Dirichlet boundary conditions
#c.constrain(2., mesh.facesLeft)
#c.constrain(1., mesh.facesRight)

## Neumann boundary conditions
c.faceGrad.constrain(-1, where=mesh.facesLeft)
c.faceGrad.constrain( -c.faceValue , where=mesh.facesRight)

Eq = DiffusionTerm(coeff=1.0)
Eq.cacheMatrix()
Eq.cacheRHSvector()
Eq.solve(var=c)
m = Eq.matrix.numpyArray
b = Eq.RHSvector

此代码无法解决，但我确实可以看到矩阵和 RHS：

m= array([[-2.,  2.,  0.],
          [ 2., -4.,  2.],
          [ 0.,  2., -2.]])

b= array([-1. ,  0. ,  0.5])

矩阵m显然是奇异的，因为源项未包含在最后一行中。关于如何包含它的任何建议？

Answer 1

[编辑：添加一阶实现的推导和演示]

一般边界条件存在已知问题。

可以实现这样的边界条件作为源。使用DiffusionTerm 的离散化， $\sum_f D_f (n\cdot\nabla(y))_f A_f$我们将$f=R$其视为特殊情况并替代所需的边界条件-n\cdot\nabla(y) - y = 0。我们通过D_(f=R)在DiffusionTerm

c.faceGrad.constrain([-1], where=mesh.facesLeft)

D = 1.
Dface = FaceVariable(mesh=mesh, value=D)
Dface.setValue(0., where=mesh.facesRight)

然后添加一个隐式源D_(f=R) (n\cdot\nabla(y))_(f=R) A_(f=R)或D_(f=R) (-y)_(f=R) A_(f=R). 源是在单元中心定义的，因此我们将单元定位在$f=R$

mask_coeff = (mesh.facesRight * mesh.faceNormals).divergence

然后添加源

Af = mesh._faceAreas[mesh.facesRight.value][0]
Eq = DiffusionTerm(coeff=Dface) - ImplicitSourceTerm(coeff=mask_coeff * D * Af)

这种处理只有 0 阶精度，因为它对ImplicitSourceTerm单元中心的值进行操作，而边界条件是在相邻的面中心定义的。

我们可以通过将单元格值沿边界条件的梯度投影到面上来使边界条件在空间中一阶准确： y_(f=R) ~ y_P + n\cdot\nabla(y)_(f=R) dPR其中y_P是y单元格中心最接近的值f=R，是点到面dPR的距离。PR

因此边界条件-n\cdot\nabla(y)_(f=R) - y_(f=R) = 0变为-n\cdot\nabla(y)_(f=R) - (y_P + n\cdot\nabla(y)_(f=R) dPR) = 0，我们可以求解n\cdot\nabla(y)_(f=R) = -y_P / (1 + dPR)。因此，对应于边界的隐式源DiffusionTerm是D_(f=R) (-y_P / (1 + dPR)) A_(f=R)

dPR = mesh._cellDistances[mesh.facesRight.value][0]
Af = mesh._faceAreas[mesh.facesRight.value][0]
Eq = DiffusionTerm(coeff=Dface) - ImplicitSourceTerm(coeff=mask_coeff * D * Af / (1 + dPR))

这是去年夏天在 FiPy 邮件列表中讨论的，从https://www.mail-archive.com/fipy@nist.gov/msg03671.html开始。是的，现在一切都很笨拙。