【问题标题】:FiPy convection with a given velocity field给定速度场的 FiPy 对流
【发布时间】:2019-11-04 14:04:30
【问题描述】:

我正在使用 FiPy 对给定速度场中的化学物质的对流/扩散进行建模,但在将此场设置为 ConvectionTerm 系数时遇到了问题。

我的速度场由二维数组表示,假设Ugrid 表示水平速度,Vgrid 表示垂直速度,表示每个单元面的值。出于这个原因,我想说将此字段设置为 ConvectionTerm 的 coeff 参数的正确方法是将其分配给 FaceVariable。

但是,我不知道如何将这两个数组作为 value 传递给 FaceVariable。显然,我通过使用value = [np.ravel(Ugrid),np.ravel(Vgrid)] 将字段设置为 CellVariable 没有问题,并且对流模拟似乎也很有意义,但我认为这不正确,正如我在上面简要提到的那样。

有什么建议吗?

【问题讨论】:

    标签: python pde fipy


    【解决方案1】:

    使用 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_interpxy 内,否则它会给出 nan 值。一旦我们有了插值,我们就可以设置 FiPy 速度场。

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

    注意ConvectionTerm 的系数可以是FaceVariableCellVariable。一旦我们有了速度,我们就可以建立和求解方程了。

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

    这对我来说运行没有错误。

    Link to complete script for this problem

    【讨论】:

      猜你喜欢
      • 1970-01-01
      • 2018-04-20
      • 2017-06-22
      • 1970-01-01
      • 1970-01-01
      • 1970-01-01
      • 1970-01-01
      • 2017-12-02
      • 1970-01-01
      相关资源
      最近更新 更多