Общие граничные условия

Буду признателен за помощь в установке общих граничных условий, -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 явно сингулярна, потому что исходный термин не был включен в последнюю строку. Любые предложения о том, как включить его?


fipy
person Scott Calabrese Barton    schedule 13.03.2017    source источник

Ответы (1)


arrow_upward
0
arrow_downward

[править: добавить вывод и демонстрацию реализации 1-го порядка]

Существуют известные проблемы с общими граничными условиями.

Такое граничное условие можно реализовать в качестве источника. Используя дискретизацию DiffusionTerm $\sum_f D_f (n\cdot\nabla(y))_f A_f$ мы рассматривать $f=R$ как частный случай и заменить желаемое граничное условие -n\cdot\nabla(y) - y = 0. Мы достигаем этого, обнуляя D_(f=R) в 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) (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 работает со значением в центре ячейки, тогда как граничное условие определяется в центре соседней грани.

Мы можем сделать граничное условие 1-го порядка точным в пространстве, спроецировав значение ячейки на грань по градиенту от граничного условия: 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. Да, это все довольно неуклюже прямо сейчас.

person jeguyer    schedule 14.03.2017
comment
Спасибо, кажется, стало ближе. Однако мне непонятно, как вводится термин потока, возможно, это как-то связано с атрибутом дивергенции? Этот подход также, по-видимому, ограничивает номинальные значения и значения центра ячейки равными на границе. (См. этот график для приведенного выше примера.) Это вводит ошибку, которая масштабируется с размером ячейки. - person Scott Calabrese Barton; 15.03.2017
comment
У меня было ощущение, что вы захотите узнать, почему источник выглядит именно так. Хороший повод для меня, чтобы выяснить, почему это выглядит так. Я отредактировал свой ответ, чтобы описать вывод и дать лучшую экстраполяцию на границы. - person jeguyer; 17.03.2017