FiPy конвекция с заданным полем скоростей - PullRequest
Новая
       73

FiPy конвекция с заданным полем скоростей

0 голосов
/ 04 ноября 2019

Я использую FiPy для моделирования конвекции / диффузии химического вещества в заданном поле скоростей, и у меня возникают проблемы с установкой этого поля в качестве коэффициента ConvectionTerm.

Мое поле скоростей выражается двумя двумерными массивамискажем, Ugrid для горизонтальных скоростей и Vgrid для вертикальных, представляющих значения на гранях каждой ячейки. По этой причине я бы сказал, что правильный способ установить это поле в качестве параметра coeff для ConvectionTerm - это присвоить его FaceVariable.

Однако я не понимаю, как передать эти два массива как value для FaceVariable. Очевидно, у меня нет проблем с установкой поля в качестве CellVariable с помощью value = [np.ravel(Ugrid),np.ravel(Vgrid)], и симуляция конвекции также имеет смысл, но я не думаю, что это правильно, как я кратко упомянул выше.

Есть предложения?

1 Ответ

1 голос
/ 13 ноября 2019

Важным соображением для FiPy является то, что индексация сетки не основана на индексации сетки, поэтому при сопоставлении массива сетки со значениями ячеек не следует предполагать, что сетка упорядочена каким-либо конкретным образом. Это делает необходимым использование некоторой формы интерполяции. Давайте сначала построим поле скорости с индексацией сетки, где точки лежат на точках сетки (не центрах ячеек). 

from fipy import Grid2D, CellVariable, FaceVariable
from fipy import ConvectionTerm, DiffusionTerm, TransientTerm
import numpy as np
from scipy import interpolate

dx = 0.5
nx = 7
dy = 0.5
ny = 7

xy = np.zeros((nx, ny, 2))
xy[:, :,  0] = np.mgrid[0:nx, 0:ny][0] * dx
xy[:, :, 1] = np.mgrid[0:nx, 0:ny][1] * dy

u = xy[..., 0] + xy[..., 1]
v = xy[..., 0] * xy[..., 1]

Теперь у нас есть поле скорости (u, v), каждое из которых имеет форму nx, ny и соответствующий емукоординаты xy, формы (nx, ny, 2). Допустим, мы хотим интерполировать это на сетку FiPy с тем же доменом, но с другой сеткой. 

m = Grid2D(nx=3, ny=3, dx=1.0, dy=1.0)

Сетка не обязательно должна совпадать с точками сетки для поля скорости. Затем мы можем интерполировать эту сетку с помощью 

xy_interp = np.array(m.faceCenters).swapaxes(0, 1)
u_interp = interpolate.griddata(xy.reshape(-1, 2), u.flatten(), xy_interp, method='cubic')
v_interp = interpolate.griddata(xy.reshape(-1, 2), v.flatten(), xy_interp, method='cubic')

, где xy_interp - это грани центра сетки. Обратите внимание, что использование griddata требует, чтобы xy_interp находилось внутри xy, иначе это дает значения nan. Получив интерполированные значения, мы можем установить поле скорости FiPy. 

velocity = FaceVariable(mesh=m, rank=1)
velocity[0, :] = u_interp
velocity[1, :] = v_interp

Обратите внимание, что коэффициент для ConvectionTerm может быть FaceVariable или CellVariable. Как только у нас будет скорость, мы готовы установить и решить уравнение. 

var = CellVariable(mesh=m)
eqn = (TransientTerm() + ConvectionTerm(velocity) == DiffusionTerm(1.0))
eqn.solve(var, dt=1.0)

Для меня это работает без ошибок.

Ссылка на полный сценарий для этой задачи

...