具有给定速度场的 FiPy 对流
FiPy convection with a given velocity field
我正在使用 FiPy 对给定速度场中的化学物质进行 convection/diffusion 建模,但我在将此场设置为 ConvectionTerm 系数时遇到问题。
我的速度场由两个二维数组表示,假设水平速度为 Ugrid
,垂直速度为 Vgrid
,表示每个单元格面上的值。出于这个原因,我认为将此字段设置为 ConvectionTerm 的 coeff
参数的正确方法是将其分配给 FaceVariable。
但是,我不知道如何将这两个数组作为 value
传递给 FaceVariable。显然,我在使用 value = [np.ravel(Ugrid),np.ravel(Vgrid)]
将字段设置为 CellVariable 时没有问题,对流模拟似乎也有意义,但我认为这是不正确的,正如我在上面简要提到的那样。
有什么建议吗?
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)
这对我来说没有错误。
我正在使用 FiPy 对给定速度场中的化学物质进行 convection/diffusion 建模,但我在将此场设置为 ConvectionTerm 系数时遇到问题。
我的速度场由两个二维数组表示,假设水平速度为 Ugrid
,垂直速度为 Vgrid
,表示每个单元格面上的值。出于这个原因,我认为将此字段设置为 ConvectionTerm 的 coeff
参数的正确方法是将其分配给 FaceVariable。
但是,我不知道如何将这两个数组作为 value
传递给 FaceVariable。显然,我在使用 value = [np.ravel(Ugrid),np.ravel(Vgrid)]
将字段设置为 CellVariable 时没有问题,对流模拟似乎也有意义,但我认为这是不正确的,正如我在上面简要提到的那样。
有什么建议吗?
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)
这对我来说没有错误。