2

如何在 fipy 网格中将法线方向的通量显式设置为特定值，而不限制面内的通量分量？

Neumann 边界条件可以指定为：(1) 垂直于边界面的通量的固定分量，或 (2) 作为面处通量的完整规范。默认的fipy条件是前者（值= 0），但显式方法（faceGrad.constrain）是后者。通过将以下代码添加到 fipy diffusion.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]
4

3 回答 3

2

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

于 2016-09-19T14:29:30.643 回答
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）设置面梯度对解决方案没有影响。

于 2016-09-19T16:49:12.690 回答
1

关于固定通量边界条件的讨论确实回答了我的问题，但我没有理解非常简短的文档。下面是一个工作示例，说明了如何在一个简单的 2D 案例中应用它，类似于 mesh20x20 示例。

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()
于 2016-09-26T23:24:40.610 回答