SUNTANS模型学习（5）——学习cavity算例

最新推荐文章于 2022-11-06 11:42:27 发布

HHU_HuLihan

最新推荐文章于 2022-11-06 11:42:27 发布

阅读量407

点赞数

分类专栏：非静压模型 SUNTANS 代码编译文章标签：学习算法

本文链接：https://blog.csdn.net/weixin_42848392/article/details/126818035

版权

非静压模型同时被 3 个专栏收录

30 篇文章 4 订阅

订阅专栏

SUNTANS

10 篇文章 0 订阅

订阅专栏

代码编译

8 篇文章 0 订阅

订阅专栏

学习cavity算例

简介
网格配置
定解条件设置
- 初始条件
- 边界条件
参数配置
模拟结果
参考文献

简介

Cavity Flow（方腔驱动流）是CFD的经典算例。该算例可用于测试SUNTANS的动量平流过程（advection of momentum），无滑移边界（no-slip）条件。
在一个正方形区域内，其三个壁面被设置为无滑移边界，剩余一个壁面做匀速平移以驱动方腔内的流动。该算例文件夹（/examples/cavity），包含了多套网格配置：

工况配置	网格类型	运行模型所用的执行命令	备注
1	三角形网格（X-Y平面）	make testXY	-
2	矩形网格（X-Y平面）	make testXY-quad	-
3	三角形网格（X-Z平面）	make testXZ	-
4	矩形网格（X-X平面）	make testXZ-quad	-

关于该算例的更多内容可参考 Zang et. al (1994) ¹:。

网格配置

网格配置如下图所示：

三角形网格（X-Y平面）
共4016个网格。
矩形网格（X-Y平面），网格数为30×30，其x和y方向的网格大小分布均满足Chebychev多项式

注：上述两个网格对应的case均是2D模拟，故我们设定垂向网格总数为Nkmax=1。

三角形网格（X-Z平面）
网格的平面分布如下图所示；垂向均分为128个等距层，初始水深为1.0m。
矩形网格（X-Z平面）
水平网格示意图和三维网格示意图如下所示；垂向分为30个不等距层，层高分布与水平网格大小分布相似，满足Chebychev多项式；初始水深为1.0m。

定解条件设置

通过Makefile我们可以得知，以下各个工况分别调用了如下的代码文件、参数文件：

工况配置	网格类型	代码文件	参数文件
1	三角形网格（X-Y平面）	boundariesXY.c, initialization-constant-dz.c	cellsXY.dat, edgesXY.dat, pointsXY.dat, suntansXY.dat
2	矩形网格（X-Y平面）	boundariesXY.c, initialization-constant-dz.c	cellsXY-quad.dat, edgesXY-quad.dat, pointsXY-quad.dat, suntansXY-quad.dat
3	三角形网格（X-Z平面）	boundariesXZ.c, initialization-constant-dz.c	cellsXZ.dat, edgesXZ.dat, pointsXZ.dat, suntansXZ.dat
4	矩形网格（X-X平面）	boundariesXY.c, initialization-chebychev-dz.c	cellsXZ-quad.dat, edgesXZ-quad.dat, pointsXZ-quad.dat, suntansXZ-quad.dat

上述代码文件对应了模型的定解条件。

初始条件

对于1、2两个工况，计算域在X-Y平面上，初始化均采用initialization-constant-dz.c。在这个c语言文件中，我们主要关注两个子函数：

/*
 * Function: GetDZ
 * Usage: GetDZ(dz,depth,Nkmax,myproc);
 * ------------------------------------
 * Returns the vertical grid spacing in the array dz.
 *
 */
int GetDZ(REAL *dz, REAL depth, REAL localdepth, int Nkmax, int myproc) {
  int k;
  REAL zk, zkp1, *z = (REAL *)SunMalloc(Nkmax*sizeof(REAL),"GetDZ");
  
  if(Nkmax==1) {
    dz[0]=depth;
  } else {
    for(k=0;k<Nkmax;k++) {
      dz[k]=1.0/Nkmax;
    }
  }
  return Nkmax;
}
  
/*
 * Function: ReturnDepth
 * Usage: grid->dv[n]=ReturnDepth(grid->xv[n],grid->yv[n]);
 * --------------------------------------------------------
 * Helper function to create a bottom bathymetry.  Used in
 * grid.c in the GetDepth function when IntDepth is 0.
 *
 */
REAL ReturnDepth(REAL x, REAL y) {
  return 1.0;
}

子函数GetDZ表示，当Nkmax=1时，模型垂向层数为1，水深为depth；depth的值则在函数ReturnDepth中得到，即depth=1.0。此外，初始的速度、盐度、温度均设置为0。

对于工况3，初始化仍采用initialization-constant-dz.c。此时Nkmax=128，子函数GetDZ显示模型的垂向网格将被均分。

对于工况4，初始化采用initialization-chebychev-dz.c。相比于initialization-constant-dz.c，此时的子函数GetDZ有着些许的不同。

int GetDZ(REAL *dz, REAL depth, REAL localdepth, int Nkmax, int myproc) {
  int k;
  REAL zk, zkp1, *z = (REAL *)SunMalloc(Nkmax*sizeof(REAL),"GetDZ");
  
  if(Nkmax==1) {
    dz[0]=depth;
  } else {
    for(k=0;k<Nkmax;k++) {
      zk=-depth*0.5*(1-cos(((REAL)k)*PI/(REAL)Nkmax));
      zkp1=-depth*0.5*(1-cos(((REAL)(k+1))*PI/(REAL)Nkmax));
      dz[k]=zk-zkp1;
    }
  }
  return Nkmax;
}

在工况4中，Nkmax=30，垂向网格每层的层高满足如下的公式：
$\Delta z_k = -0.5 h_0 (1-cos(k \pi \times Nkmax)) + 0.5 h_0 (1-cos((k+1) \pi \times Nkmax)), k =0,1,2 ... Nkmax-1$
此外，初始的速度、盐度、温度也均设置为0。

边界条件

对于工况1、2，计算域在X-Y平面上，边界条件通过boundariesXY.c指定。在这个c语言文件中，我们重点关注以下子函数：

/*
 * Function: BoundaryVelocities
 * Usage: BoundaryVelocities(grid,phys,prop,myproc);
 * -------------------------------------------------
 * This will set the values of u,v,w, and h at the boundaries.
 * 
 */
void BoundaryVelocities(gridT *grid, physT *phys, propT *prop, int myproc, MPI_Comm comm) {
  int j, k, jptr, ib, boundary_index;
  REAL tau=0.5;

  for(jptr=grid->edgedist[4];jptr<grid->edgedist[5];jptr++) {
    j = grid->edgep[jptr];
    ib=grid->grad[2*j];
    boundary_index = jptr-grid->edgedist[2];

    for(k=grid->ctop[ib];k<grid->Nk[ib];k++) {
      if(grid->ye[j]<0.000010) 
        phys->boundary_u[boundary_index][k]= 1.0;
      else
       phys->boundary_u[boundary_index][k] = 0.0;
       phys->boundary_v[boundary_index][k] = 0.0;
       phys->boundary_w[boundary_index][k] = 0.5*(phys->w[ib][k]+phys->w[ib][k+1]);
    }
  }
}

在输入文件中，所有边界的类型均被设为了mark=4，它们可以通过edgedist[4]~edgedist[5]来检索。对于y=0的南边界，边界水平速度boundary_u为1.0m/s，其余边界的速度均为0。

对于工况3、4，计算域在X-Z平面上，边界条件通过boundariesXZ.c指定。其子函数BoundaryVelocities内容如下：

/*
 * Function: BoundaryVelocities
 * Usage: BoundaryVelocities(grid,phys,prop,myproc);
 * -------------------------------------------------
 * This will set the values of u,v,w, and h at the boundaries.
 * 
 */
void BoundaryVelocities(gridT *grid, physT *phys, propT *prop, int myproc, MPI_Comm comm) {
  int i, j, k, iptr, jptr, ib, boundary_index;
  REAL tau=0.5;

  for(jptr=grid->edgedist[4];jptr<grid->edgedist[5];jptr++) {
    
    j = grid->edgep[jptr];
    ib=grid->grad[2*j];
    boundary_index = jptr-grid->edgedist[2];

    for(k=grid->ctop[ib];k<grid->Nk[ib];k++) {
      if(grid->xe[j]<0.5) 
        phys->boundary_w[boundary_index][k]= -1.0;
      phys->boundary_v[boundary_index][k] = 0.0;
      phys->boundary_u[boundary_index][k] = 0;
    }
  }
}

在此，移动边界被设置在x=0，即模型的西边界；其垂向速度boundary_w被设为了-1.0m/s（方向沿z方向向下）。

参数配置

模型参数通过suntans.dat中，如上表所示，对于不同的工况，有四种不同的suntans.dat文件。一些重要的参数设置如下：

对流项的离散均采用中心差分格式（nonlinear = 2）；
在XY工况种采用静压模拟（nonhydrostatic = 0），在XZ工况中采用非静压模拟（nonhydrostatic = 1）；
流体粘度根据Re=3200来设置（即nu = 0.0003125、nu_H = 0.0003125）；

模拟结果

以下展示运行至流场稳定时，三角形网格（X-Y平面）算例所对应的流速大小云图与流速矢量图。
以下展示运行至流场稳定时，矩形网格（X-Y平面）算例所对应的流速大小云图与流速矢量图。
以下展示运行至流场稳定时，三角形网格（X-Z平面）算例所对应的流速大小云图与流速矢量图。
以下展示运行至流场稳定时，矩形网格（X-Z平面）算例所对应的流速大小云图与流速矢量图。

参考文献

Zang Y , Street R L , Koseff J R . A Non-staggered Grid, Fractional Step Method for Time-Dependent Incompressible Navier-Stokes Equations in Curvilinear Coordinates[J]. J.comp.physics, 1994, 114(1):18-33. ↩︎