一般边界条件

Question

我希望能帮助我设置一般边界条件，-grad(y) + g(y) = 0 其中g 是未知y 的某个函数。这是一个我无法开始工作的简单 1D 示例：

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

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

有known issues 具有一般边界条件。

可以实现这样的边界条件作为源。使用discretization of the DiffusionTerm$\sum_f D_f (n\cdot\nabla(y))_f A_f$，我们将$f=R$ 视为一种特殊情况并替代所需的边界条件-n\cdot\nabla(y) - y = 0。我们通过在DiffusionTerm 中将D_(f=R) 归零来实现这一点

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 是最靠近f=R 的单元格中心处y 的值，dPR 是点P 到面R 的距离。

因此边界条件-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 开始。是的，现在一切都很笨拙。