2016-09-18 221 views
2

我如何明確地設置fipy網格中的邊界面的流量法線爲一個特定值,而不會限制面內的流量分量? (1)通量的固定分量垂直於邊界面,或(2)作爲面的通量的完整規定。默認的fipy條件是前者(值= 0),但顯式方法(faceGrad.constrain)是後者。通過將以下代碼添加到fipy diffusion.mesh20x20示例的結尾並注意不同的臉部漸變結果可以理解此問題。您如何在Fipy中指定Neumann(法向面的固定通量)邊界條件?

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] 

回答

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)對解決方案沒有影響。

+0

如果您希望一個邊界條件取決於另一個邊界的梯度的大小,並且僅僅需要設置正常分量,該怎麼辦?說,'phi1.constrain(fp.numerix.dot(phi2.faceGrad,phi2.faceGrad))'我們只想在邊界上指定phi2的正常分量。這個答案是否意味着這樣做不起作用? (我沒有一個很好的例子,只是好奇)。也許@ jeguyer的回答在這種情況下是適當的? – muon

+1

@ muon;我認爲約束方法對於大多數任意複雜的相互依賴的邊界條件都是失敗的,並且需要源項方法(@jeguyer答案)。 – wd15

+0

有道理。謝謝! – muon

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()