如何在Fipy中指定诺伊曼(法向面固定通量)边界条件

Question

如何显式地将五边形网格中边界面的通量法线设置为特定值，而不约束面内的通量分量?

诺依曼边界条件可以表示为:(1)法向边界面通量的固定分量，或(2)边界面上通量的完整规格。默认的固定条件是前者(value = 0)，但显式方法(facegrad . constraint)是后者。可以通过将以下代码添加到fipy扩散的末尾来理解这个问题。Mesh20x20的例子，并注意不同的人脸梯度结果。

facesNeumann = mesh.exteriorFaces & ~facesTopLeft & ~facesBottomRight
print 'grad(phi) with default Neumann BC...'
print phi.faceGrad.value.T[facesNeumann.value]
phi.faceGrad.constrain(0.,facesNeumann)
DiffusionTerm().solve(var=phi)
print 'and with explicit Neumann BC...'
print phi.faceGrad.value.T[facesNeumann.value]

Answer 1

请参见关于固定通量边界条件的讨论。基本上，您可以将包含所需边界通量散度的源添加到FiPy的默认无通量条件中。

Answer 2

对于独立于求解方程的变量，仅仅指定法向边界的通量，似乎没有任何方法可以做到这一点。例如，语法可能是phi.faceGrad[0].constrain(...)，但目前在FiPy中不起作用。对于任意定向的面，也可能很难实现。

然而，出于实际目的，在FiPy中求解方程时不使用与边界相切的分量，只有法向分量对解有任何影响。例如，使用以下代码

import fipy as fp
mesh = fp.Grid2D(nx=2, ny=1)
var = fp.CellVariable(mesh=mesh)
var.constrain(1, mesh.facesLeft)
var.constrain(0, mesh.facesRight)
#var.faceGrad.constrain(0, mesh.facesTop)
fp.DiffusionTerm().solve(var)
print 'face gradient on top plane:',var.faceGrad[0, mesh.facesTop.value]
print 'variable value:',var

输出

face gradient on top plane: [-0.5 -0.5]
variable value: [ 0.75  0.25]

答案是正确的，但是顶部的面梯度是-0.5。但是，当#var.faceGrad.constrain(0, mesh.facesTop)行没有注释时，结果是

face gradient on top plane: [ 0.  0.]
variable value: [ 0.75  0.25]

切向面梯度现在是0，但答案是一样的。关键是切向设置面梯度(通过.constrain)对溶液没有影响。

Answer 3

关于固定通量边界条件的讨论确实回答了我的问题，但是我没有理解文档，它非常简短。下面是一个有效的例子，说明了如何在一个简单的2D情况下应用它，类似于网格20x20的例子。

import matplotlib.pyplot as plt
from fipy import *
nx = 20
ny = nx
dx = 1.
dy = dx
L = dx * nx
msh = Grid2D(dx=dx, dy=dy, nx=nx, ny=ny)
xFace, yFace = msh.faceCenters
xCell,yCell = msh.cellCenters    
phi = CellVariable(name = "solution variable",
                   mesh = msh,
                   value = 0.)    
D = FaceVariable(name='diffusion coefficient',mesh=msh,value=1.)
# Dirichlet condition on top face
valueTop = FaceVariable(name='valueTop',mesh=msh,value=xFace*0.1-1)
phi.constrain(valueTop, msh.facesTop)
# Fixed flux (Neumann) on base
D.constrain(0., msh.facesBottom)
fluxBottom = -0.05
eq = DiffusionTerm(coeff=D) + (msh.facesBottom * fluxBottom).divergence
eq.solve(var=phi)
# confirm only y component of gradient is zero
print phi.faceGrad.value.T[msh.facesBottom.value]
plt.scatter(phi.value, yCell)
plt.show()
viewer = Viewer(vars=phi, datamin=-1., datamax=1.)
viewer.plot()