一般边界条件

Question

我希望能帮助你设置一般的边界条件，其中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))

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