一般邊界條件

我很感謝設置一般邊界條件的幫助，-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，顯然是單數的，因爲源詞不包含在最後一行。有關如何包含它的任何建議？

來源

2017-03-13 Scott Calabrese Barton

[編輯：添加推導和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

c.faceGrad.constrain([-1], where=mesh.facesLeft) 

D = 1. 
Dface = FaceVariable(mesh=mesh, value=D) 
Dface.setValue(0., where=mesh.facesRight)

歸零D_(f=R)，然後加上一個隱含的源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是y在最接近小區中心到f=R和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/[email protected]/msg03671.html。是的，現在都很笨拙。

來源

2017-03-14 17:02:33 jeguyer

謝謝，這似乎越來越接近。然而，我不清楚流量術語是如何引入的，也許它與分歧屬性有關係？這種方法似乎也限制了邊界處的面值和小區中心值相等。（對於上面的例子，請參見[this plot]（http://imgur.com/a/jJz3M））。這引入了一個隨單元大小而變化的錯誤。 –

我有一種感覺，你會想知道爲什麼源代碼看起來像那樣。對我來說，找出它爲什麼是這樣的好藉口。我編輯了我的答案來描述推導，並給出了更好的邊界推斷。 – jeguyer

一般邊界條件

回答

相關問題